﻿using System;
using System.Collections.Generic;
using System.Linq;
using SpiceNet.Circuits;
using SpiceNet.Diagnostics;
using SpiceNet.Models.BSIM3p;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    /// <summary>
    /// An instance of the BSIM3 Mosfet
    /// </summary>
    public class BSIM3 : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, IP> ParamTable => _pTable;

        #region Parameters
        private enum _c
        {
            /* device parameters */
            BSIM3_W = 1,
            BSIM3_L = 2,
            BSIM3_AS = 3,
            BSIM3_AD = 4,
            BSIM3_PS = 5,
            BSIM3_PD = 6,
            BSIM3_NRS = 7,
            BSIM3_NRD = 8,
            BSIM3_OFF = 9,
            BSIM3_IC_VBS = 10,
            BSIM3_IC_VDS = 11,
            BSIM3_IC_VGS = 12,
            BSIM3_IC = 13,
            BSIM3_NQSMOD = 14,
            BSIM3_ACNQSMOD = 15,
            /* device questions */
            BSIM3_DNODE = 751,
            BSIM3_GNODE = 752,
            BSIM3_SNODE = 753,
            BSIM3_BNODE = 754,
            BSIM3_DNODEPRIME = 755,
            BSIM3_SNODEPRIME = 756,
            BSIM3_VBD = 757,
            BSIM3_VBS = 758,
            BSIM3_VGS = 759,
            BSIM3_VDS = 760,
            BSIM3_CD = 761,
            BSIM3_CBS = 762,
            BSIM3_CBD = 763,
            BSIM3_GM = 764,
            BSIM3_GDS = 765,
            BSIM3_GMBS = 766,
            BSIM3_GBD = 767,
            BSIM3_GBS = 768,
            BSIM3_QB = 769,
            BSIM3_CQB = 770,
            BSIM3_QG = 771,
            BSIM3_CQG = 772,
            BSIM3_QD = 773,
            BSIM3_CQD = 774,
            BSIM3_CGG = 775,
            BSIM3_CGD = 776,
            BSIM3_CGS = 777,
            BSIM3_CBG = 778,
            BSIM3_CAPBD = 779,
            BSIM3_CQBD = 780,
            BSIM3_CAPBS = 781,
            BSIM3_CQBS = 782,
            BSIM3_CDG = 783,
            BSIM3_CDD = 784,
            BSIM3_CDS = 785,
            BSIM3_VON = 786,
            BSIM3_VDSAT = 787,
            BSIM3_QBS = 788,
            BSIM3_QBD = 789,
            BSIM3_SOURCECONDUCT = 790,
            BSIM3_DRAINCONDUCT = 791,
            BSIM3_CBDB = 792,
            BSIM3_CBSB = 793
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "l", new IP(IF.IOP, (int)_c.BSIM3_L, T.REAL, "Length") },
            { "w", new IP(IF.IOP, (int)_c.BSIM3_W, T.REAL, "Width") },
            { "ad", new IP(IF.IOP, (int)_c.BSIM3_AD, T.REAL, "Drain area") },
            { "as", new IP(IF.IOP, (int)_c.BSIM3_AS, T.REAL, "Source area") },
            { "pd", new IP(IF.IOP, (int)_c.BSIM3_PD, T.REAL, "Drain perimeter") },
            { "ps", new IP(IF.IOP, (int)_c.BSIM3_PS, T.REAL, "Source perimeter") },
            { "nrd", new IP(IF.IOP, (int)_c.BSIM3_NRD, T.REAL, "Number of squares in drain") },
            { "nrs", new IP(IF.IOP, (int)_c.BSIM3_NRS, T.REAL, "Number of squares in source") },
            { "off", new IP(IF.IOP, (int)_c.BSIM3_OFF, T.FLAG, "Device is initially off") },
            { "nqsmod", new IP(IF.IOP, (int)_c.BSIM3_NQSMOD, T.INTEGER, "Non-quasi-static model selector") },
            { "acnqsmod", new IP(IF.IOP, (int)_c.BSIM3_ACNQSMOD, T.INTEGER, "AC NQS model selector") },
            { "ic", new IP(IF.IP, (int)_c.BSIM3_IC, T.REALVEC, "Vector of DS,GS,BS initial voltages") },
            { "gmbs", new IP(IF.OP, (int)_c.BSIM3_GMBS, T.REAL, "Gmb") },
            { "gm", new IP(IF.OP, (int)_c.BSIM3_GM, T.REAL, "Gm") },
            { "gds", new IP(IF.OP, (int)_c.BSIM3_GDS, T.REAL, "Gds") },
            { "vdsat", new IP(IF.OP, (int)_c.BSIM3_VDSAT, T.REAL, "Vdsat") },
            { "vth", new IP(IF.OP, (int)_c.BSIM3_VON, T.REAL, "Vth") },
            { "id", new IP(IF.OP, (int)_c.BSIM3_CD, T.REAL, "Ids") },
            { "vbs", new IP(IF.OP, (int)_c.BSIM3_VBS, T.REAL, "Vbs") },
            { "vgs", new IP(IF.OP, (int)_c.BSIM3_VGS, T.REAL, "Vgs") },
            { "vds", new IP(IF.OP, (int)_c.BSIM3_VDS, T.REAL, "Vds") }
        };

        public Parameter<double> BSIM3drainArea { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3drainPerimeter { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3drainSquares { get; } = new Parameter<double>(1.0);
        public Parameter<double> BSIM3icVBS { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3icVDS { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3icVGS { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3l { get; } = new Parameter<double>(5.0e-6);
        public Parameter<double> BSIM3sourceArea { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3sourcePerimeter { get; } = new Parameter<double>(0.0);
        public Parameter<double> BSIM3sourceSquares { get; } = new Parameter<double>(1.0);
        public Parameter<double> BSIM3w { get; } = new Parameter<double>(5.0e-6);
        public Parameter<int> BSIM3nqsMod { get; } = new Parameter<int>(0);
        public Parameter<int> BSIM3acnqsMod { get; } = new Parameter<int>(0);
        public bool BSIM3off { get; set; } = false;
        #endregion

        #region Constants
        private const int BSIM3vbd = 0;
        private const int BSIM3vbs = 1;
        private const int BSIM3vgs = 2;
        private const int BSIM3vds = 3;

        private const int BSIM3qb = 4;
        private const int BSIM3cqb = 5;
        private const int BSIM3qg = 6;
        private const int BSIM3cqg = 7;
        private const int BSIM3qd = 8;
        private const int BSIM3cqd = 9;

        private const int BSIM3qbs = 10;
        private const int BSIM3qbd = 11;

        private const int BSIM3qcheq = 12;
        private const int BSIM3cqcheq = 13;
        private const int BSIM3qcdump = 14;
        private const int BSIM3cqcdump = 15;

        private const int BSIM3qdef = 16;

        private const double DELTA_1 = 0.02;
        private const double DELTA_2 = 0.02;
        private const double DELTA_3 = 0.02;
        private const double DELTA_4 = 0.02;
        #endregion

        /// <summary>
        /// Private variables for this instance
        /// </summary>
        private int BSIM3dNode, BSIM3gNode, BSIM3sNode, BSIM3bNode, BSIM3dNodePrime, BSIM3sNodePrime, BSIM3qNode;
        private int BSIM3state;

        //  MCJ 
        private double BSIM3ueff, BSIM3thetavth, BSIM3von, BSIM3vdsat, BSIM3cgdo, BSIM3cgso, BSIM3vjsm, BSIM3IsEvjsm, BSIM3vjdm, BSIM3IsEvjdm,
            BSIM3sourceConductance, BSIM3drainConductance;
        private int BSIM3mode;

        //  OP point 
        private double BSIM3qinv, BSIM3cd, BSIM3cbs, BSIM3cbd, BSIM3csub, BSIM3gm, BSIM3gds, BSIM3gmbs, BSIM3gbd, BSIM3gbs, BSIM3gbbs, BSIM3gbgs,
            BSIM3gbds, BSIM3cggb, BSIM3cgdb, BSIM3cgsb, BSIM3cbgb, BSIM3cbdb, BSIM3cbsb, BSIM3cdgb, BSIM3cddb, BSIM3cdsb, BSIM3capbd, BSIM3capbs,
            BSIM3cqgb, BSIM3cqdb, BSIM3cqsb, BSIM3cqbb, BSIM3qgate, BSIM3qbulk, BSIM3qdrn, BSIM3gtau, BSIM3gtg, BSIM3gtd, BSIM3gts, BSIM3gtb,
            BSIM3rds,  //  Noise bugfix 
            BSIM3Vgsteff, BSIM3Vdseff, BSIM3Abulk, BSIM3AbovVgst2Vtm, BSIM3taunet;
        private SizeDependParam param = null;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public BSIM3(string name) : base(name, 4, typeof(BSIM3Model)) { }

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        /// <param name="drainnode"></param>
        /// <param name="gatenode"></param>
        /// <param name="sourcenode"></param>
        /// <param name="bulknode"></param>
        public BSIM3(string name, string drainnode, string gatenode, string sourcenode, string bulknode) : this(name)
        {
            Connect(drainnode, gatenode, sourcenode, bulknode);
        }

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        /// <param name="dnode">The drain node</param>
        /// <param name="gnode">The gate node</param>
        /// <param name="snode">The source node</param>
        /// <param name="bnode">The bulk node</param>
        /// <param name="w">The width of the transistor</param>
        /// <param name="l">The length of the transistor</param>
        public BSIM3(string name, string dnode, string gnode, string snode, string bnode, double w = 5.0e-6, double l = 5.0e-6) : this(name)
        {
            Connect(dnode, gnode, snode, bnode);
            BSIM3w.Par(w);
            BSIM3l.Par(l);
        }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            double[] v;
            switch ((_c)id)
            {
                case _c.BSIM3_W:
                    BSIM3w.Par((double)value);
                    break;
                case _c.BSIM3_L:
                    BSIM3l.Par((double)value);
                    break;
                case _c.BSIM3_AS:
                    BSIM3sourceArea.Par((double)value);
                    break;
                case _c.BSIM3_AD:
                    BSIM3drainArea.Par((double)value);
                    break;
                case _c.BSIM3_PS:
                    BSIM3sourcePerimeter.Par((double)value);
                    break;
                case _c.BSIM3_PD:
                    BSIM3drainPerimeter.Par((double)value);
                    break;
                case _c.BSIM3_NRS:
                    BSIM3sourceSquares.Par((double)value);
                    break;
                case _c.BSIM3_NRD:
                    BSIM3drainSquares.Par((double)value);
                    break;
                case _c.BSIM3_OFF:
                    BSIM3off = (bool)value;
                    break;
                case _c.BSIM3_IC_VBS:
                    BSIM3icVBS.Par((double)value);
                    break;
                case _c.BSIM3_IC_VDS:
                    BSIM3icVDS.Par((double)value);
                    break;
                case _c.BSIM3_IC_VGS:
                    BSIM3icVGS.Par((double)value);
                    break;
                case _c.BSIM3_NQSMOD:
                    BSIM3nqsMod.Par((int)value);
                    break;
                case _c.BSIM3_ACNQSMOD:
                    BSIM3acnqsMod.Par((int)value);
                    break;
                case _c.BSIM3_IC:
                    v = (double[])value;
                    switch (v.Length)
                    {
                        case 3:
                            BSIM3icVBS.Par(v[2]); goto case 2;
                        case 2:
                            BSIM3icVGS.Par(v[1]); goto case 1;
                        case 1:
                            BSIM3icVDS.Par(v[0]); break;
                        default:
                            throw new BadParameterException();
                    }
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Get parameter
        /// </summary>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.BSIM3_L:
                    return BSIM3l.Value;
                case _c.BSIM3_W:
                    return BSIM3w.Value;
                case _c.BSIM3_AS:
                    return BSIM3sourceArea.Value;
                case _c.BSIM3_AD:
                    return BSIM3drainArea.Value;
                case _c.BSIM3_PS:
                    return BSIM3sourcePerimeter.Value;
                case _c.BSIM3_PD:
                    return BSIM3drainPerimeter.Value;
                case _c.BSIM3_NRS:
                    return BSIM3sourceSquares.Value;
                case _c.BSIM3_NRD:
                    return BSIM3drainSquares.Value;
                case _c.BSIM3_OFF:
                    return BSIM3off;
                case _c.BSIM3_NQSMOD:
                    return BSIM3nqsMod.Value;
                case _c.BSIM3_ACNQSMOD:
                    return BSIM3acnqsMod.Value;
                case _c.BSIM3_IC_VBS:
                    return BSIM3icVBS.Value;
                case _c.BSIM3_IC_VDS:
                    return BSIM3icVDS.Value;
                case _c.BSIM3_IC_VGS:
                    return BSIM3icVGS.Value;
                case _c.BSIM3_DNODE:
                    return BSIM3dNode;
                case _c.BSIM3_GNODE:
                    return BSIM3gNode;
                case _c.BSIM3_SNODE:
                    return BSIM3sNode;
                case _c.BSIM3_BNODE:
                    return BSIM3bNode;
                case _c.BSIM3_DNODEPRIME:
                    return BSIM3dNodePrime;
                case _c.BSIM3_SNODEPRIME:
                    return BSIM3sNodePrime;
                case _c.BSIM3_SOURCECONDUCT:
                    return BSIM3sourceConductance;
                case _c.BSIM3_DRAINCONDUCT:
                    return BSIM3drainConductance;
                case _c.BSIM3_VBD:
                    return ckt.State.States[0][BSIM3state + BSIM3vbd];
                case _c.BSIM3_VBS:
                    return ckt.State.States[0][BSIM3state + BSIM3vbs];
                case _c.BSIM3_VGS:
                    return ckt.State.States[0][BSIM3state + BSIM3vgs];
                case _c.BSIM3_VDS:
                    return ckt.State.States[0][BSIM3state + BSIM3vds];
                case _c.BSIM3_CD:
                    return BSIM3cd;
                case _c.BSIM3_CBS:
                    return BSIM3cbs;
                case _c.BSIM3_CBD:
                    return BSIM3cbd;
                case _c.BSIM3_GM:
                    return BSIM3gm;
                case _c.BSIM3_GDS:
                    return BSIM3gds;
                case _c.BSIM3_GMBS:
                    return BSIM3gmbs;
                case _c.BSIM3_GBD:
                    return BSIM3gbd;
                case _c.BSIM3_GBS:
                    return BSIM3gbs;
                case _c.BSIM3_QB:
                    return ckt.State.States[0][BSIM3state + BSIM3qb];
                case _c.BSIM3_CQB:
                    return ckt.State.States[0][BSIM3state + BSIM3cqb];
                case _c.BSIM3_QG:
                    return ckt.State.States[0][BSIM3state + BSIM3qg];
                case _c.BSIM3_CQG:
                    return ckt.State.States[0][BSIM3state + BSIM3cqg];
                case _c.BSIM3_QD:
                    return ckt.State.States[0][BSIM3state + BSIM3qd];
                case _c.BSIM3_CQD:
                    return ckt.State.States[0][BSIM3state + BSIM3cqd];
                case _c.BSIM3_CGG:
                    return BSIM3cggb;
                case _c.BSIM3_CGD:
                    return BSIM3cgdb;
                case _c.BSIM3_CGS:
                    return BSIM3cgsb;
                case _c.BSIM3_CDG:
                    return BSIM3cdgb;
                case _c.BSIM3_CDD:
                    return BSIM3cddb;
                case _c.BSIM3_CDS:
                    return BSIM3cdsb;
                case _c.BSIM3_CBG:
                    return BSIM3cbgb;
                case _c.BSIM3_CBDB:
                    return BSIM3cbdb;
                case _c.BSIM3_CBSB:
                    return BSIM3cbsb;
                case _c.BSIM3_CAPBD:
                    return BSIM3capbd;
                case _c.BSIM3_CAPBS:
                    return BSIM3capbs;
                case _c.BSIM3_VON:
                    return BSIM3von;
                case _c.BSIM3_VDSAT:
                    return BSIM3vdsat;
                case _c.BSIM3_QBS:
                    return ckt.State.States[0][BSIM3state + BSIM3qbs];
                case _c.BSIM3_QBD:
                    return ckt.State.States[0][BSIM3state + BSIM3qbd];
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Setup the BSIM3 Mosfet
        /// </summary>
        /// <param name="ckt"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            BSIM3Model model = mod as BSIM3Model;

            // Handle NQS
            if (!BSIM3acnqsMod.Given)
                BSIM3acnqsMod.Value = model.BSIM3acnqsMod;
            else if ((BSIM3acnqsMod != 0) && (BSIM3acnqsMod != 1))
            {
                BSIM3acnqsMod.Value = model.BSIM3acnqsMod;
                CircuitWarning.Warning(this, string.Format("Warning: acnqsMod has been set to its global value {0}", model.BSIM3acnqsMod.Value));
            }

            // Add extra nodes if the bipolar transistor has conductances at the collector or emitter
            List<NodeType> extra = new List<NodeType>();

            // NQS mode - add charge
            bool addCharge = false;
            if (BSIM3nqsMod != 0)
            {
                extra.Add(NodeType.Voltage);
                addCharge = true;
            }

            bool addSourcePrime = false;
            if ((model.BSIM3sheetResistance > 0.0) && (BSIM3drainSquares > 0.0))
            {
                extra.Add(NodeType.Voltage);
                addSourcePrime = true;
            }

            bool addDrainPrime = false;
            if ((model.BSIM3sheetResistance > 0.0) && (BSIM3sourceSquares > 0.0))
            {
                extra.Add(NodeType.Voltage);
                addDrainPrime = true;
            }

            // Bind the nodes
            CircuitNode[] nodes = BindNodes(ckt, extra.ToArray());
            int index = 4;
            BSIM3dNode = nodes[0].Number;
            BSIM3gNode = nodes[1].Number;
            BSIM3sNode = nodes[2].Number;
            BSIM3bNode = nodes[3].Number;
            BSIM3qNode = addCharge ? nodes[index++].Number : 0;
            BSIM3dNodePrime = addDrainPrime ? nodes[index++].Number : BSIM3dNode;
            BSIM3sNodePrime = addSourcePrime ? nodes[index++].Number : BSIM3sNode;

            // Reserve states
            BSIM3state = CreateState(ckt, 17);

            // TODO: Sensitivity
        }

        /// <summary>
        /// Temperature-dependent calculations
        /// </summary>
        /// <param name="ckt"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            BSIM3Model model = mod as BSIM3Model;
            BSIM3Model.TempParams tp = p as BSIM3Model.TempParams;
            double tmp;
            double tmp1;
            double tmp2;
            double tmp3;
            double T2;
            double T3;
            double T4;
            double T5;
            double Ldrn;
            double Wdrn;
            double Inv_L;
            double Inv_W;
            double Inv_LW;
            double Nvtm;
            double SourceSatCurrent;
            double DrainSatCurrent;

            // Create a new size entry if necessary
            Tuple<double, double> size_param = new Tuple<double, double>(BSIM3w.Value, BSIM3l.Value);
            if (model.Params.ContainsKey(size_param))
                this.param = model.Params[size_param];
            else
            {
                this.param = new SizeDependParam();
                model.Params.Add(size_param, this.param);

                Ldrn = BSIM3l;
                Wdrn = BSIM3w;
                param.Length = Ldrn;
                param.Width = Wdrn;

                tp.T0 = Math.Pow(Ldrn, model.BSIM3Lln);
                tp.T1 = Math.Pow(Wdrn, model.BSIM3Lwn);
                tmp1 = model.BSIM3Ll / tp.T0 + model.BSIM3Lw / tp.T1
                       + model.BSIM3Lwl / (tp.T0 * tp.T1);
                param.BSIM3dl = model.BSIM3Lint + tmp1;
                tmp2 = model.BSIM3Llc / tp.T0 + model.BSIM3Lwc / tp.T1
                       + model.BSIM3Lwlc / (tp.T0 * tp.T1);
                param.BSIM3dlc = model.BSIM3dlc + tmp2;

                T2 = Math.Pow(Ldrn, model.BSIM3Wln);
                T3 = Math.Pow(Wdrn, model.BSIM3Wwn);
                tmp1 = model.BSIM3Wl / T2 + model.BSIM3Ww / T3
                       + model.BSIM3Wwl / (T2 * T3);
                param.BSIM3dw = model.BSIM3Wint + tmp1;
                tmp2 = model.BSIM3Wlc / T2 + model.BSIM3Wwc / T3
                       + model.BSIM3Wwlc / (T2 * T3);
                param.BSIM3dwc = model.BSIM3dwc + tmp2;

                param.BSIM3leff = BSIM3l - 2.0 * param.BSIM3dl;
                if (param.BSIM3leff <= 0.0)
                    throw new CircuitException(String.Format("BSIM3: mosfet {0}, model {1}: Effective channel length <= 0", Name, model.Name));

                param.BSIM3weff = BSIM3w - 2.0 * param.BSIM3dw;
                if (param.BSIM3weff <= 0.0)
                    throw new CircuitException(String.Format("BSIM3: mosfet {0}, model {1}: Effective channel width <= 0", Name, model.Name));

                param.BSIM3leffCV = BSIM3l - 2.0 * param.BSIM3dlc;
                if (param.BSIM3leffCV <= 0.0)
                    throw new CircuitException(String.Format("BSIM3: mosfet {0}, model {1}: Effective channel length for C-V <= 0", Name, model.Name));

                param.BSIM3weffCV = BSIM3w - 2.0 * param.BSIM3dwc;
                if (param.BSIM3weffCV <= 0.0)
                    throw new CircuitException(String.Format("BSIM3: mosfet {0}, model {1}: Effective channel width for C-V <= 0", Name, model.Name));

                if (model.BSIM3binUnit == 1)
                {
                    Inv_L = 1.0e-6 / param.BSIM3leff;
                    Inv_W = 1.0e-6 / param.BSIM3weff;
                    Inv_LW = 1.0e-12 / (param.BSIM3leff
                                        * param.BSIM3weff);
                }
                else
                {
                    Inv_L = 1.0 / param.BSIM3leff;
                    Inv_W = 1.0 / param.BSIM3weff;
                    Inv_LW = 1.0 / (param.BSIM3leff
                                    * param.BSIM3weff);
                }
                param.BSIM3cdsc = model.BSIM3cdsc
                                    + model.BSIM3lcdsc * Inv_L
                                    + model.BSIM3wcdsc * Inv_W
                                    + model.BSIM3pcdsc * Inv_LW;
                param.BSIM3cdscb = model.BSIM3cdscb
                                     + model.BSIM3lcdscb * Inv_L
                                     + model.BSIM3wcdscb * Inv_W
                                     + model.BSIM3pcdscb * Inv_LW;

                param.BSIM3cdscd = model.BSIM3cdscd
                                     + model.BSIM3lcdscd * Inv_L
                                     + model.BSIM3wcdscd * Inv_W
                                     + model.BSIM3pcdscd * Inv_LW;

                param.BSIM3cit = model.BSIM3cit
                                   + model.BSIM3lcit * Inv_L
                                   + model.BSIM3wcit * Inv_W
                                   + model.BSIM3pcit * Inv_LW;
                param.BSIM3nfactor = model.BSIM3nfactor
                                       + model.BSIM3lnfactor * Inv_L
                                       + model.BSIM3wnfactor * Inv_W
                                       + model.BSIM3pnfactor * Inv_LW;
                param.BSIM3xj = model.BSIM3xj
                                  + model.BSIM3lxj * Inv_L
                                  + model.BSIM3wxj * Inv_W
                                  + model.BSIM3pxj * Inv_LW;
                param.BSIM3vsat = model.BSIM3vsat
                                    + model.BSIM3lvsat * Inv_L
                                    + model.BSIM3wvsat * Inv_W
                                    + model.BSIM3pvsat * Inv_LW;
                param.BSIM3at = model.BSIM3at
                                  + model.BSIM3lat * Inv_L
                                  + model.BSIM3wat * Inv_W
                                  + model.BSIM3pat * Inv_LW;
                param.BSIM3a0 = model.BSIM3a0
                                  + model.BSIM3la0 * Inv_L
                                  + model.BSIM3wa0 * Inv_W
                                  + model.BSIM3pa0 * Inv_LW;

                param.BSIM3ags = model.BSIM3ags
                                   + model.BSIM3lags * Inv_L
                                   + model.BSIM3wags * Inv_W
                                   + model.BSIM3pags * Inv_LW;

                param.BSIM3a1 = model.BSIM3a1
                                  + model.BSIM3la1 * Inv_L
                                  + model.BSIM3wa1 * Inv_W
                                  + model.BSIM3pa1 * Inv_LW;
                param.BSIM3a2 = model.BSIM3a2
                                  + model.BSIM3la2 * Inv_L
                                  + model.BSIM3wa2 * Inv_W
                                  + model.BSIM3pa2 * Inv_LW;
                param.BSIM3keta = model.BSIM3keta
                                    + model.BSIM3lketa * Inv_L
                                    + model.BSIM3wketa * Inv_W
                                    + model.BSIM3pketa * Inv_LW;
                param.BSIM3nsub = model.BSIM3nsub
                                    + model.BSIM3lnsub * Inv_L
                                    + model.BSIM3wnsub * Inv_W
                                    + model.BSIM3pnsub * Inv_LW;
                param.BSIM3npeak = model.BSIM3npeak
                                     + model.BSIM3lnpeak * Inv_L
                                     + model.BSIM3wnpeak * Inv_W
                                     + model.BSIM3pnpeak * Inv_LW;
                param.BSIM3ngate = model.BSIM3ngate
                                     + model.BSIM3lngate * Inv_L
                                     + model.BSIM3wngate * Inv_W
                                     + model.BSIM3pngate * Inv_LW;
                param.BSIM3gamma1 = model.BSIM3gamma1
                                      + model.BSIM3lgamma1 * Inv_L
                                      + model.BSIM3wgamma1 * Inv_W
                                      + model.BSIM3pgamma1 * Inv_LW;
                param.BSIM3gamma2 = model.BSIM3gamma2
                                      + model.BSIM3lgamma2 * Inv_L
                                      + model.BSIM3wgamma2 * Inv_W
                                      + model.BSIM3pgamma2 * Inv_LW;
                param.BSIM3vbx = model.BSIM3vbx
                                   + model.BSIM3lvbx * Inv_L
                                   + model.BSIM3wvbx * Inv_W
                                   + model.BSIM3pvbx * Inv_LW;
                param.BSIM3vbm = model.BSIM3vbm
                                   + model.BSIM3lvbm * Inv_L
                                   + model.BSIM3wvbm * Inv_W
                                   + model.BSIM3pvbm * Inv_LW;
                param.BSIM3xt = model.BSIM3xt
                                  + model.BSIM3lxt * Inv_L
                                  + model.BSIM3wxt * Inv_W
                                  + model.BSIM3pxt * Inv_LW;
                param.BSIM3vfb = model.BSIM3vfb
                                   + model.BSIM3lvfb * Inv_L
                                   + model.BSIM3wvfb * Inv_W
                                   + model.BSIM3pvfb * Inv_LW;
                param.BSIM3k1 = model.BSIM3k1
                                  + model.BSIM3lk1 * Inv_L
                                  + model.BSIM3wk1 * Inv_W
                                  + model.BSIM3pk1 * Inv_LW;
                param.BSIM3kt1 = model.BSIM3kt1
                                   + model.BSIM3lkt1 * Inv_L
                                   + model.BSIM3wkt1 * Inv_W
                                   + model.BSIM3pkt1 * Inv_LW;
                param.BSIM3kt1l = model.BSIM3kt1l
                                    + model.BSIM3lkt1l * Inv_L
                                    + model.BSIM3wkt1l * Inv_W
                                    + model.BSIM3pkt1l * Inv_LW;
                param.BSIM3k2 = model.BSIM3k2
                                  + model.BSIM3lk2 * Inv_L
                                  + model.BSIM3wk2 * Inv_W
                                  + model.BSIM3pk2 * Inv_LW;
                param.BSIM3kt2 = model.BSIM3kt2
                                   + model.BSIM3lkt2 * Inv_L
                                   + model.BSIM3wkt2 * Inv_W
                                   + model.BSIM3pkt2 * Inv_LW;
                param.BSIM3k3 = model.BSIM3k3
                                  + model.BSIM3lk3 * Inv_L
                                  + model.BSIM3wk3 * Inv_W
                                  + model.BSIM3pk3 * Inv_LW;
                param.BSIM3k3b = model.BSIM3k3b
                                   + model.BSIM3lk3b * Inv_L
                                   + model.BSIM3wk3b * Inv_W
                                   + model.BSIM3pk3b * Inv_LW;
                param.BSIM3w0 = model.BSIM3w0
                                  + model.BSIM3lw0 * Inv_L
                                  + model.BSIM3ww0 * Inv_W
                                  + model.BSIM3pw0 * Inv_LW;
                param.BSIM3nlx = model.BSIM3nlx
                                   + model.BSIM3lnlx * Inv_L
                                   + model.BSIM3wnlx * Inv_W
                                   + model.BSIM3pnlx * Inv_LW;
                param.BSIM3dvt0 = model.BSIM3dvt0
                                    + model.BSIM3ldvt0 * Inv_L
                                    + model.BSIM3wdvt0 * Inv_W
                                    + model.BSIM3pdvt0 * Inv_LW;
                param.BSIM3dvt1 = model.BSIM3dvt1
                                    + model.BSIM3ldvt1 * Inv_L
                                    + model.BSIM3wdvt1 * Inv_W
                                    + model.BSIM3pdvt1 * Inv_LW;
                param.BSIM3dvt2 = model.BSIM3dvt2
                                    + model.BSIM3ldvt2 * Inv_L
                                    + model.BSIM3wdvt2 * Inv_W
                                    + model.BSIM3pdvt2 * Inv_LW;
                param.BSIM3dvt0w = model.BSIM3dvt0w
                                     + model.BSIM3ldvt0w * Inv_L
                                     + model.BSIM3wdvt0w * Inv_W
                                     + model.BSIM3pdvt0w * Inv_LW;
                param.BSIM3dvt1w = model.BSIM3dvt1w
                                     + model.BSIM3ldvt1w * Inv_L
                                     + model.BSIM3wdvt1w * Inv_W
                                     + model.BSIM3pdvt1w * Inv_LW;
                param.BSIM3dvt2w = model.BSIM3dvt2w
                                     + model.BSIM3ldvt2w * Inv_L
                                     + model.BSIM3wdvt2w * Inv_W
                                     + model.BSIM3pdvt2w * Inv_LW;
                param.BSIM3drout = model.BSIM3drout
                                     + model.BSIM3ldrout * Inv_L
                                     + model.BSIM3wdrout * Inv_W
                                     + model.BSIM3pdrout * Inv_LW;
                param.BSIM3dsub = model.BSIM3dsub
                                    + model.BSIM3ldsub * Inv_L
                                    + model.BSIM3wdsub * Inv_W
                                    + model.BSIM3pdsub * Inv_LW;
                param.BSIM3vth0 = model.BSIM3vth0
                                    + model.BSIM3lvth0 * Inv_L
                                    + model.BSIM3wvth0 * Inv_W
                                    + model.BSIM3pvth0 * Inv_LW;
                param.BSIM3ua = model.BSIM3ua
                                  + model.BSIM3lua * Inv_L
                                  + model.BSIM3wua * Inv_W
                                  + model.BSIM3pua * Inv_LW;
                param.BSIM3ua1 = model.BSIM3ua1
                                   + model.BSIM3lua1 * Inv_L
                                   + model.BSIM3wua1 * Inv_W
                                   + model.BSIM3pua1 * Inv_LW;
                param.BSIM3ub = model.BSIM3ub
                                  + model.BSIM3lub * Inv_L
                                  + model.BSIM3wub * Inv_W
                                  + model.BSIM3pub * Inv_LW;
                param.BSIM3ub1 = model.BSIM3ub1
                                   + model.BSIM3lub1 * Inv_L
                                   + model.BSIM3wub1 * Inv_W
                                   + model.BSIM3pub1 * Inv_LW;
                param.BSIM3uc = model.BSIM3uc
                                  + model.BSIM3luc * Inv_L
                                  + model.BSIM3wuc * Inv_W
                                  + model.BSIM3puc * Inv_LW;
                param.BSIM3uc1 = model.BSIM3uc1
                                   + model.BSIM3luc1 * Inv_L
                                   + model.BSIM3wuc1 * Inv_W
                                   + model.BSIM3puc1 * Inv_LW;
                param.BSIM3u0 = model.BSIM3u0
                                  + model.BSIM3lu0 * Inv_L
                                  + model.BSIM3wu0 * Inv_W
                                  + model.BSIM3pu0 * Inv_LW;
                param.BSIM3ute = model.BSIM3ute
                                   + model.BSIM3lute * Inv_L
                                   + model.BSIM3wute * Inv_W
                                   + model.BSIM3pute * Inv_LW;
                param.BSIM3voff = model.BSIM3voff
                                    + model.BSIM3lvoff * Inv_L
                                    + model.BSIM3wvoff * Inv_W
                                    + model.BSIM3pvoff * Inv_LW;
                param.BSIM3delta = model.BSIM3delta
                                     + model.BSIM3ldelta * Inv_L
                                     + model.BSIM3wdelta * Inv_W
                                     + model.BSIM3pdelta * Inv_LW;
                param.BSIM3rdsw = model.BSIM3rdsw
                                    + model.BSIM3lrdsw * Inv_L
                                    + model.BSIM3wrdsw * Inv_W
                                    + model.BSIM3prdsw * Inv_LW;
                param.BSIM3prwg = model.BSIM3prwg
                                    + model.BSIM3lprwg * Inv_L
                                    + model.BSIM3wprwg * Inv_W
                                    + model.BSIM3pprwg * Inv_LW;
                param.BSIM3prwb = model.BSIM3prwb
                                    + model.BSIM3lprwb * Inv_L
                                    + model.BSIM3wprwb * Inv_W
                                    + model.BSIM3pprwb * Inv_LW;
                param.BSIM3prt = model.BSIM3prt
                                   + model.BSIM3lprt * Inv_L
                                   + model.BSIM3wprt * Inv_W
                                   + model.BSIM3pprt * Inv_LW;
                param.BSIM3eta0 = model.BSIM3eta0
                                    + model.BSIM3leta0 * Inv_L
                                    + model.BSIM3weta0 * Inv_W
                                    + model.BSIM3peta0 * Inv_LW;
                param.BSIM3etab = model.BSIM3etab
                                    + model.BSIM3letab * Inv_L
                                    + model.BSIM3wetab * Inv_W
                                    + model.BSIM3petab * Inv_LW;
                param.BSIM3pclm = model.BSIM3pclm
                                    + model.BSIM3lpclm * Inv_L
                                    + model.BSIM3wpclm * Inv_W
                                    + model.BSIM3ppclm * Inv_LW;
                param.BSIM3pdibl1 = model.BSIM3pdibl1
                                      + model.BSIM3lpdibl1 * Inv_L
                                      + model.BSIM3wpdibl1 * Inv_W
                                      + model.BSIM3ppdibl1 * Inv_LW;
                param.BSIM3pdibl2 = model.BSIM3pdibl2
                                      + model.BSIM3lpdibl2 * Inv_L
                                      + model.BSIM3wpdibl2 * Inv_W
                                      + model.BSIM3ppdibl2 * Inv_LW;
                param.BSIM3pdiblb = model.BSIM3pdiblb
                                      + model.BSIM3lpdiblb * Inv_L
                                      + model.BSIM3wpdiblb * Inv_W
                                      + model.BSIM3ppdiblb * Inv_LW;
                param.BSIM3pscbe1 = model.BSIM3pscbe1
                                      + model.BSIM3lpscbe1 * Inv_L
                                      + model.BSIM3wpscbe1 * Inv_W
                                      + model.BSIM3ppscbe1 * Inv_LW;
                param.BSIM3pscbe2 = model.BSIM3pscbe2
                                      + model.BSIM3lpscbe2 * Inv_L
                                      + model.BSIM3wpscbe2 * Inv_W
                                      + model.BSIM3ppscbe2 * Inv_LW;
                param.BSIM3pvag = model.BSIM3pvag
                                    + model.BSIM3lpvag * Inv_L
                                    + model.BSIM3wpvag * Inv_W
                                    + model.BSIM3ppvag * Inv_LW;
                param.BSIM3wr = model.BSIM3wr
                                  + model.BSIM3lwr * Inv_L
                                  + model.BSIM3wwr * Inv_W
                                  + model.BSIM3pwr * Inv_LW;
                param.BSIM3dwg = model.BSIM3dwg
                                   + model.BSIM3ldwg * Inv_L
                                   + model.BSIM3wdwg * Inv_W
                                   + model.BSIM3pdwg * Inv_LW;
                param.BSIM3dwb = model.BSIM3dwb
                                   + model.BSIM3ldwb * Inv_L
                                   + model.BSIM3wdwb * Inv_W
                                   + model.BSIM3pdwb * Inv_LW;
                param.BSIM3b0 = model.BSIM3b0
                                  + model.BSIM3lb0 * Inv_L
                                  + model.BSIM3wb0 * Inv_W
                                  + model.BSIM3pb0 * Inv_LW;
                param.BSIM3b1 = model.BSIM3b1
                                  + model.BSIM3lb1 * Inv_L
                                  + model.BSIM3wb1 * Inv_W
                                  + model.BSIM3pb1 * Inv_LW;
                param.BSIM3alpha0 = model.BSIM3alpha0
                                      + model.BSIM3lalpha0 * Inv_L
                                      + model.BSIM3walpha0 * Inv_W
                                      + model.BSIM3palpha0 * Inv_LW;
                param.BSIM3alpha1 = model.BSIM3alpha1
                                      + model.BSIM3lalpha1 * Inv_L
                                      + model.BSIM3walpha1 * Inv_W
                                      + model.BSIM3palpha1 * Inv_LW;
                param.BSIM3beta0 = model.BSIM3beta0
                                     + model.BSIM3lbeta0 * Inv_L
                                     + model.BSIM3wbeta0 * Inv_W
                                     + model.BSIM3pbeta0 * Inv_LW;
                /* CV model */
                param.BSIM3elm = model.BSIM3elm
                                   + model.BSIM3lelm * Inv_L
                                   + model.BSIM3welm * Inv_W
                                   + model.BSIM3pelm * Inv_LW;
                param.BSIM3cgsl = model.BSIM3cgsl
                                    + model.BSIM3lcgsl * Inv_L
                                    + model.BSIM3wcgsl * Inv_W
                                    + model.BSIM3pcgsl * Inv_LW;
                param.BSIM3cgdl = model.BSIM3cgdl
                                    + model.BSIM3lcgdl * Inv_L
                                    + model.BSIM3wcgdl * Inv_W
                                    + model.BSIM3pcgdl * Inv_LW;
                param.BSIM3ckappa = model.BSIM3ckappa
                                      + model.BSIM3lckappa * Inv_L
                                      + model.BSIM3wckappa * Inv_W
                                      + model.BSIM3pckappa * Inv_LW;
                param.BSIM3cf = model.BSIM3cf
                                  + model.BSIM3lcf * Inv_L
                                  + model.BSIM3wcf * Inv_W
                                  + model.BSIM3pcf * Inv_LW;
                param.BSIM3clc = model.BSIM3clc
                                   + model.BSIM3lclc * Inv_L
                                   + model.BSIM3wclc * Inv_W
                                   + model.BSIM3pclc * Inv_LW;
                param.BSIM3cle = model.BSIM3cle
                                   + model.BSIM3lcle * Inv_L
                                   + model.BSIM3wcle * Inv_W
                                   + model.BSIM3pcle * Inv_LW;
                param.BSIM3vfbcv = model.BSIM3vfbcv
                                     + model.BSIM3lvfbcv * Inv_L
                                     + model.BSIM3wvfbcv * Inv_W
                                     + model.BSIM3pvfbcv * Inv_LW;
                param.BSIM3acde = model.BSIM3acde
                                    + model.BSIM3lacde * Inv_L
                                    + model.BSIM3wacde * Inv_W
                                    + model.BSIM3pacde * Inv_LW;
                param.BSIM3moin = model.BSIM3moin
                                    + model.BSIM3lmoin * Inv_L
                                    + model.BSIM3wmoin * Inv_W
                                    + model.BSIM3pmoin * Inv_LW;
                param.BSIM3noff = model.BSIM3noff
                                    + model.BSIM3lnoff * Inv_L
                                    + model.BSIM3wnoff * Inv_W
                                    + model.BSIM3pnoff * Inv_LW;
                param.BSIM3voffcv = model.BSIM3voffcv
                                      + model.BSIM3lvoffcv * Inv_L
                                      + model.BSIM3wvoffcv * Inv_W
                                      + model.BSIM3pvoffcv * Inv_LW;

                param.BSIM3abulkCVfactor = 1.0 + Math.Pow((param.BSIM3clc
                                                        / param.BSIM3leffCV),
                                                       param.BSIM3cle);

                tp.T0 = (tp.TRatio - 1.0);
                param.BSIM3ua = param.BSIM3ua + param.BSIM3ua1 * tp.T0;
                param.BSIM3ub = param.BSIM3ub + param.BSIM3ub1 * tp.T0;
                param.BSIM3uc = param.BSIM3uc + param.BSIM3uc1 * tp.T0;
                if (param.BSIM3u0 > 1.0)
                    param.BSIM3u0 = param.BSIM3u0 / 1.0e4;

                param.BSIM3u0temp = param.BSIM3u0
                                      * Math.Pow(tp.TRatio, param.BSIM3ute);
                param.BSIM3vsattemp = param.BSIM3vsat - param.BSIM3at
                                        * tp.T0;
                param.BSIM3rds0 = (param.BSIM3rdsw + param.BSIM3prt * tp.T0)
                                    / Math.Pow(param.BSIM3weff * 1E6, param.BSIM3wr);

                BSIM3checkModel(ckt, model);

                param.BSIM3cgdo = (model.BSIM3cgdo + param.BSIM3cf)
                                    * param.BSIM3weffCV;
                param.BSIM3cgso = (model.BSIM3cgso + param.BSIM3cf)
                                    * param.BSIM3weffCV;
                param.BSIM3cgbo = model.BSIM3cgbo * param.BSIM3leffCV;

                tp.T0 = param.BSIM3leffCV * param.BSIM3leffCV;
                param.BSIM3tconst = param.BSIM3u0temp * param.BSIM3elm / (model.BSIM3cox
                                      * param.BSIM3weffCV * param.BSIM3leffCV * tp.T0);

                if (!model.BSIM3npeak.Given && model.BSIM3gamma1.Given)
                {
                    tp.T0 = param.BSIM3gamma1 * model.BSIM3cox;
                    param.BSIM3npeak = 3.021E22 * tp.T0 * tp.T0;
                }

                param.BSIM3phi = 2.0 * tp.Vtm0
                                   * Math.Log(param.BSIM3npeak / tp.ni);

                param.BSIM3sqrtPhi = Math.Sqrt(param.BSIM3phi);
                param.BSIM3phis3 = param.BSIM3sqrtPhi * param.BSIM3phi;

                param.BSIM3Xdep0 = Math.Sqrt(2.0 * BSIM3Model.EPSSI / (BSIM3Model.Charge_q
                                          * param.BSIM3npeak * 1.0e6))
                                     * param.BSIM3sqrtPhi;
                param.BSIM3sqrtXdep0 = Math.Sqrt(param.BSIM3Xdep0);
                param.BSIM3litl = Math.Sqrt(3.0 * param.BSIM3xj
                                         * model.BSIM3tox);
                param.BSIM3vbi = tp.Vtm0 * Math.Log(1.0e20
                                              * param.BSIM3npeak / (tp.ni * tp.ni));
                param.BSIM3cdep0 = Math.Sqrt(BSIM3Model.Charge_q * BSIM3Model.EPSSI
                                          * param.BSIM3npeak * 1.0e6 / 2.0
                                          / param.BSIM3phi);

                param.BSIM3ldeb = Math.Sqrt(BSIM3Model.EPSSI * tp.Vtm0 / (BSIM3Model.Charge_q
                                         * param.BSIM3npeak * 1.0e6)) / 3.0;
                param.BSIM3acde *= Math.Pow((param.BSIM3npeak / 2.0e16), -0.25);


                if (model.BSIM3k1.Given || model.BSIM3k2.Given)
                {
                    if (!model.BSIM3k1.Given)
                    {
                        CircuitWarning.Warning(this, "Warning: k1 should be specified with k2.\n");
                        param.BSIM3k1 = 0.53;
                    }
                    if (!model.BSIM3k2.Given)
                    {
                        CircuitWarning.Warning(this, "Warning: k2 should be specified with k1.\n");
                        param.BSIM3k2 = -0.0186;
                    }
                    if (model.BSIM3nsub.Given)
                        CircuitWarning.Warning(this, "Warning: nsub is ignored because k1 or k2 is given.\n");
                    if (model.BSIM3xt.Given)
                        CircuitWarning.Warning(this, "Warning: xt is ignored because k1 or k2 is given.\n");
                    if (model.BSIM3vbx.Given)
                        CircuitWarning.Warning(this, "Warning: vbx is ignored because k1 or k2 is given.\n");
                    if (model.BSIM3gamma1.Given)
                        CircuitWarning.Warning(this, "Warning: gamma1 is ignored because k1 or k2 is given.\n");
                    if (model.BSIM3gamma2.Given)
                        CircuitWarning.Warning(this, "Warning: gamma2 is ignored because k1 or k2 is given.\n");
                }
                else
                {
                    if (!model.BSIM3vbx.Given)
                        param.BSIM3vbx = param.BSIM3phi - 7.7348e-4
                                           * param.BSIM3npeak
                                           * param.BSIM3xt * param.BSIM3xt;
                    if (param.BSIM3vbx > 0.0)
                        param.BSIM3vbx = -param.BSIM3vbx;
                    if (param.BSIM3vbm > 0.0)
                        param.BSIM3vbm = -param.BSIM3vbm;

                    if (!model.BSIM3gamma1.Given)
                        param.BSIM3gamma1 = 5.753e-12
                                              * Math.Sqrt(param.BSIM3npeak)
                                              / model.BSIM3cox;
                    if (!model.BSIM3gamma2.Given)
                        param.BSIM3gamma2 = 5.753e-12
                                              * Math.Sqrt(param.BSIM3nsub)
                                              / model.BSIM3cox;

                    tp.T0 = param.BSIM3gamma1 - param.BSIM3gamma2;
                    tp.T1 = Math.Sqrt(param.BSIM3phi - param.BSIM3vbx)
                         - param.BSIM3sqrtPhi;
                    T2 = Math.Sqrt(param.BSIM3phi * (param.BSIM3phi
                                                  - param.BSIM3vbm)) - param.BSIM3phi;
                    param.BSIM3k2 = tp.T0 * tp.T1 / (2.0 * T2 + param.BSIM3vbm);
                    param.BSIM3k1 = param.BSIM3gamma2 - 2.0
                                      * param.BSIM3k2 * Math.Sqrt(param.BSIM3phi
                                                               - param.BSIM3vbm);
                }

                if (param.BSIM3k2 < 0.0)
                {
                    tp.T0 = 0.5 * param.BSIM3k1 / param.BSIM3k2;
                    param.BSIM3vbsc = 0.9 * (param.BSIM3phi - tp.T0 * tp.T0);
                    if (param.BSIM3vbsc > -3.0)
                        param.BSIM3vbsc = -3.0;
                    else if (param.BSIM3vbsc < -30.0)
                        param.BSIM3vbsc = -30.0;
                }
                else
                {
                    param.BSIM3vbsc = -30.0;
                }
                if (param.BSIM3vbsc > param.BSIM3vbm)
                    param.BSIM3vbsc = param.BSIM3vbm;

                if (!model.BSIM3vfb.Given)
                {
                    if (model.BSIM3vth0.Given)
                    {
                        param.BSIM3vfb = model.BSIM3type * param.BSIM3vth0
                                           - param.BSIM3phi - param.BSIM3k1
                                           * param.BSIM3sqrtPhi;
                    }
                    else
                    {
                        param.BSIM3vfb = -1.0;
                    }
                }
                if (!model.BSIM3vth0.Given)
                {
                    param.BSIM3vth0 = model.BSIM3type * (param.BSIM3vfb
                                                            + param.BSIM3phi + param.BSIM3k1
                                                            * param.BSIM3sqrtPhi);
                }

                param.BSIM3k1ox = param.BSIM3k1 * model.BSIM3tox
                                    / model.BSIM3toxm;
                param.BSIM3k2ox = param.BSIM3k2 * model.BSIM3tox
                                    / model.BSIM3toxm;

                tp.T1 = Math.Sqrt(BSIM3Model.EPSSI / BSIM3Model.EPSOX * model.BSIM3tox
                          * param.BSIM3Xdep0);
                tp.T0 = Math.Exp(-0.5 * param.BSIM3dsub * param.BSIM3leff / tp.T1);
                param.BSIM3theta0vb0 = (tp.T0 + 2.0 * tp.T0 * tp.T0);

                tp.T0 = Math.Exp(-0.5 * param.BSIM3drout * param.BSIM3leff / tp.T1);
                T2 = (tp.T0 + 2.0 * tp.T0 * tp.T0);
                param.BSIM3thetaRout = param.BSIM3pdibl1 * T2
                                         + param.BSIM3pdibl2;

                tmp = Math.Sqrt(param.BSIM3Xdep0);
                tmp1 = param.BSIM3vbi - param.BSIM3phi;
                tmp2 = model.BSIM3factor1 * tmp;

                tp.T0 = -0.5 * param.BSIM3dvt1w * param.BSIM3weff
                     * param.BSIM3leff / tmp2;
                if (tp.T0 > -BSIM3Model.EXP_THRESHOLD)
                {
                    tp.T1 = Math.Exp(tp.T0);
                    T2 = tp.T1 * (1.0 + 2.0 * tp.T1);
                }
                else
                {
                    tp.T1 = BSIM3Model.MIN_EXP;
                    T2 = tp.T1 * (1.0 + 2.0 * tp.T1);
                }
                tp.T0 = param.BSIM3dvt0w * T2;
                T2 = tp.T0 * tmp1;

                tp.T0 = -0.5 * param.BSIM3dvt1 * param.BSIM3leff / tmp2;
                if (tp.T0 > -BSIM3Model.EXP_THRESHOLD)
                {
                    tp.T1 = Math.Exp(tp.T0);
                    T3 = tp.T1 * (1.0 + 2.0 * tp.T1);
                }
                else
                {
                    tp.T1 = BSIM3Model.MIN_EXP;
                    T3 = tp.T1 * (1.0 + 2.0 * tp.T1);
                }
                T3 = param.BSIM3dvt0 * T3 * tmp1;

                T4 = model.BSIM3tox * param.BSIM3phi
                     / (param.BSIM3weff + param.BSIM3w0);

                tp.T0 = Math.Sqrt(1.0 + param.BSIM3nlx / param.BSIM3leff);
                T5 = param.BSIM3k1ox * (tp.T0 - 1.0) * param.BSIM3sqrtPhi
                     + (param.BSIM3kt1 + param.BSIM3kt1l / param.BSIM3leff)
                     * (tp.TRatio - 1.0);

                tmp3 = model.BSIM3type * param.BSIM3vth0
                       - T2 - T3 + param.BSIM3k3 * T4 + T5;
                param.BSIM3vfbzb = tmp3 - param.BSIM3phi - param.BSIM3k1
                                     * param.BSIM3sqrtPhi;
                /* End of vfbzb */
            }

            /* process source/drain series resistance */
            BSIM3drainConductance = model.BSIM3sheetResistance
                                          * BSIM3drainSquares;
            if (BSIM3drainConductance > 0.0)
                BSIM3drainConductance = 1.0
                                              / BSIM3drainConductance;
            else
                BSIM3drainConductance = 0.0;

            BSIM3sourceConductance = model.BSIM3sheetResistance
                                           * BSIM3sourceSquares;
            if (BSIM3sourceConductance > 0.0)
                BSIM3sourceConductance = 1.0
                                               / BSIM3sourceConductance;
            else
                BSIM3sourceConductance = 0.0;
            BSIM3cgso = param.BSIM3cgso;
            BSIM3cgdo = param.BSIM3cgdo;

            Nvtm = model.BSIM3vtm * model.BSIM3jctEmissionCoeff;
            if ((BSIM3sourceArea <= 0.0) &&
                    (BSIM3sourcePerimeter <= 0.0))
            {
                SourceSatCurrent = 1.0e-14;
            }
            else
            {
                SourceSatCurrent = BSIM3sourceArea
                                   * model.BSIM3jctTempSatCurDensity
                                   + BSIM3sourcePerimeter
                                   * model.BSIM3jctSidewallTempSatCurDensity;
            }
            if ((SourceSatCurrent > 0.0) && (model.BSIM3ijth > 0.0))
            {
                BSIM3vjsm = Nvtm * Math.Log(model.BSIM3ijth
                                             / SourceSatCurrent + 1.0);
                BSIM3IsEvjsm = SourceSatCurrent * Math.Exp(BSIM3vjsm
                                     / Nvtm);
            }

            if ((BSIM3drainArea <= 0.0) &&
                    (BSIM3drainPerimeter <= 0.0))
            {
                DrainSatCurrent = 1.0e-14;
            }
            else
            {
                DrainSatCurrent = BSIM3drainArea
                                  * model.BSIM3jctTempSatCurDensity
                                  + BSIM3drainPerimeter
                                  * model.BSIM3jctSidewallTempSatCurDensity;
            }
            if ((DrainSatCurrent > 0.0) && (model.BSIM3ijth > 0.0))
            {
                BSIM3vjdm = Nvtm * Math.Log(model.BSIM3ijth
                                             / DrainSatCurrent + 1.0);
                BSIM3IsEvjdm = DrainSatCurrent * Math.Exp(BSIM3vjdm
                                     / Nvtm);
            }
        }

        /// <summary>
        /// Check the model parameters
        /// </summary>
        /// <param name="ckt"></param>
        private void BSIM3checkModel(Circuit ckt, BSIM3Model model)
        {
            if (model.BSIM3version != "3.3.0")
                CircuitWarning.Warning(this, "Warning: This model is BSIM3.3.0; you specified a wrong version number.\n");

            if (param.BSIM3nlx < -param.BSIM3leff)
                throw new CircuitException(String.Format("Fatal: Nlx = {0} is less than -Leff.\n", param.BSIM3nlx));
            if (model.BSIM3tox <= 0.0)
                throw new CircuitException(String.Format("Fatal: Tox = {0} is not positive.\n", model.BSIM3tox));
            if (model.BSIM3toxm <= 0.0)
                throw new CircuitException(String.Format("Fatal: Toxm = {0} is not positive.\n", model.BSIM3toxm));
            if (model.BSIM3lintnoi > param.BSIM3leff / 2)
                throw new CircuitException(String.Format("Fatal: Lintnoi = {0} is too large - Leff for noise is negative.\n", model.BSIM3lintnoi));
            if (param.BSIM3npeak <= 0.0)
                throw new CircuitException(String.Format("Fatal: Nch = {0} is not positive.\n", param.BSIM3npeak));
            if (param.BSIM3nsub <= 0.0)
                throw new CircuitException(String.Format("Fatal: Nsub = {0} is not positive.\n", param.BSIM3nsub));
            if (param.BSIM3ngate < 0.0)
                throw new CircuitException(String.Format("Fatal: Ngate = {0} Ngate is not positive.\n", param.BSIM3ngate));
            if (param.BSIM3ngate > 1.0e25)
                throw new CircuitException(String.Format("Fatal: Ngate = {0} Ngate is too high\n", param.BSIM3ngate));
            if (param.BSIM3xj <= 0.0)
                throw new CircuitException(String.Format("Fatal: Xj = {0} is not positive.\n", param.BSIM3xj));
            if (param.BSIM3dvt1 < 0.0)
                throw new CircuitException(String.Format("Fatal: Dvt1 = {0} is negative.\n", param.BSIM3dvt1));
            if (param.BSIM3dvt1w < 0.0)
                throw new CircuitException(String.Format("Fatal: Dvt1w = {0} is negative.\n", param.BSIM3dvt1w));
            if (param.BSIM3w0 == -param.BSIM3weff)
                throw new CircuitException(String.Format("Fatal: (W0 + Weff) = 0 causing divided-by-zero.\n"));
            if (param.BSIM3dsub < 0.0)
                throw new CircuitException(String.Format("Fatal: Dsub = {0} is negative.\n", param.BSIM3dsub));
            if (param.BSIM3b1 == -param.BSIM3weff)
                throw new CircuitException(String.Format("Fatal: (B1 + Weff) = 0 causing divided-by-zero.\n"));
            if (param.BSIM3u0temp <= 0.0)
                throw new CircuitException(String.Format("Fatal: u0 at current temperature = {0} is not positive.\n", param.BSIM3u0temp));
            //  Check delta parameter 
            if (param.BSIM3delta < 0.0)
                throw new CircuitException(String.Format("Fatal: Delta = {0} is less than zero.\n", param.BSIM3delta));
            if (param.BSIM3vsattemp <= 0.0)
                throw new CircuitException(String.Format("Fatal: Vsat at current temperature = {0} is not positive.\n", param.BSIM3vsattemp));
            //  Check Rout parameters 
            if (param.BSIM3pclm <= 0.0)
                throw new CircuitException(String.Format("Fatal: Pclm = {0} is not positive.\n", param.BSIM3pclm));
            if (param.BSIM3drout < 0.0)
                throw new CircuitException(String.Format("Fatal: Drout = {0} is negative.\n", param.BSIM3drout));
            if (param.BSIM3pscbe2 <= 0.0)
                CircuitWarning.Warning(this, string.Format("Warning: Pscbe2 = {0} is not positive.\n", param.BSIM3pscbe2));
            if (model.BSIM3unitLengthSidewallJctCap > 0.0 || model.BSIM3unitLengthGateSidewallJctCap > 0.0)
            {
                if (BSIM3drainPerimeter < param.BSIM3weff)
                    CircuitWarning.Warning(this, string.Format("Warning: Pd = {0} is less than W.\n", BSIM3drainPerimeter.Value));
                if (BSIM3sourcePerimeter < param.BSIM3weff)
                    CircuitWarning.Warning(this, string.Format("Warning: Ps = {0} is less than W.\n", BSIM3sourcePerimeter.Value));
            }
            if (param.BSIM3noff < 0.1)
                CircuitWarning.Warning(this, string.Format("Warning: Noff = {0} is too small.\n", param.BSIM3noff));
            if (param.BSIM3noff > 4.0)
                CircuitWarning.Warning(this, string.Format("Warning: Noff = {0} is too large.\n", param.BSIM3noff));
            if (param.BSIM3voffcv < -0.5)
                CircuitWarning.Warning(this, string.Format("Warning: Voffcv = {0} is too small.\n", param.BSIM3voffcv));
            if (param.BSIM3voffcv > 0.5)
                CircuitWarning.Warning(this, string.Format("Warning: Voffcv = {0} is too large.\n", param.BSIM3voffcv));
            if (model.BSIM3ijth < 0.0)
                throw new CircuitException(String.Format("Fatal: Ijth = {0} cannot be negative.\n", model.BSIM3ijth));
            //  Check capacitance parameters 
            if (param.BSIM3clc < 0.0)
                throw new CircuitException(String.Format("Fatal: Clc = {0} is negative.\n", param.BSIM3clc));
            if (param.BSIM3moin < 5.0)
                CircuitWarning.Warning(this, string.Format("Warning: Moin = {0} is too small.\n", param.BSIM3moin));
            if (param.BSIM3moin > 25.0)
                CircuitWarning.Warning(this, string.Format("Warning: Moin = {0} is too large.\n", param.BSIM3moin));
            if (model.BSIM3capMod == 3)
            {
                if (param.BSIM3acde < 0.4)
                    CircuitWarning.Warning(this, string.Format("Warning: Acde = {0} is too small.\n", param.BSIM3acde));
                if (param.BSIM3acde > 1.6)
                    CircuitWarning.Warning(this, string.Format("Warning: Acde = {0} is too large.\n", param.BSIM3acde));
            }
            if (model.BSIM3paramChk == 1)
            {
                //  Check L and W parameters 
                if (param.BSIM3leff <= 5.0e-8)
                    CircuitWarning.Warning(this, string.Format("Warning: Leff = {0} may be too small.\n", param.BSIM3leff));
                if (param.BSIM3leffCV <= 5.0e-8)
                    CircuitWarning.Warning(this, string.Format("Warning: Leff for CV = {0} may be too small.\n", param.BSIM3leffCV));
                if (param.BSIM3weff <= 1.0e-7)
                    CircuitWarning.Warning(this, string.Format("Warning: Weff = {0} may be too small.\n", param.BSIM3weff));
                if (param.BSIM3weffCV <= 1.0e-7)
                    CircuitWarning.Warning(this, string.Format("Warning: Weff for CV = {0} may be too small.\n", param.BSIM3weffCV));
                //  Check threshold voltage parameters 
                if (param.BSIM3nlx < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Nlx = {0} is negative.\n", param.BSIM3nlx));
                if (model.BSIM3tox < 1.0e-9)
                    CircuitWarning.Warning(this, string.Format("Warning: Tox = {0} is less than 10A.\n", model.BSIM3tox));
                if (param.BSIM3npeak <= 1.0e15)
                    CircuitWarning.Warning(this, string.Format("Warning: Nch = {0} may be too small.\n", param.BSIM3npeak));
                else if (param.BSIM3npeak >= 1.0e21)
                    CircuitWarning.Warning(this, string.Format("Warning: Nch = {0} may be too large.\n", param.BSIM3npeak));
                if (param.BSIM3nsub <= 1.0e14)
                    CircuitWarning.Warning(this, string.Format("Warning: Nsub = {0} may be too small.\n", param.BSIM3nsub));
                else if (param.BSIM3nsub >= 1.0e21)
                    CircuitWarning.Warning(this, string.Format("Warning: Nsub = {0} may be too large.\n", param.BSIM3nsub));
                if ((param.BSIM3ngate > 0.0) && (param.BSIM3ngate <= 1.0e18))
                    CircuitWarning.Warning(this, string.Format("Warning: Ngate = {0} is less than 1.E18cm^-3.\n", param.BSIM3ngate));
                if (param.BSIM3dvt0 < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Dvt0 = {0} is negative.\n", param.BSIM3dvt0));
                if (Math.Abs(1.0e-6 / (param.BSIM3w0 + param.BSIM3weff)) > 10.0)
                    CircuitWarning.Warning(this, string.Format("Warning: (W0 + Weff) may be too small.\n"));
                //  Check subthreshold parameters 
                if (param.BSIM3nfactor < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Nfactor = {0} is negative.\n", param.BSIM3nfactor));
                if (param.BSIM3cdsc < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Cdsc = {0} is negative.\n", param.BSIM3cdsc));
                if (param.BSIM3cdscd < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Cdscd = {0} is negative.\n", param.BSIM3cdscd));
                //  Check DIBL parameters 
                if (param.BSIM3eta0 < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Eta0 = {0} is negative.\n", param.BSIM3eta0));
                //  Check Abulk parameters 
                if (Math.Abs(1.0e-6 / (param.BSIM3b1 + param.BSIM3weff)) > 10.0)
                    CircuitWarning.Warning(this, string.Format("Warning: (B1 + Weff) may be too small.\n"));
                //  Check Saturation parameters 
                if (param.BSIM3a2 < 0.01)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: A2 = {0} is too small. Set to 0.01.\n", param.BSIM3a2));
                    param.BSIM3a2 = 0.01;
                }
                else if (param.BSIM3a2 > 1.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: A2 = {0} is larger than 1. A2 is set to 1 and A1 is set to 0.\n", param.BSIM3a2));
                    param.BSIM3a2 = 1.0;
                    param.BSIM3a1 = 0.0;
                }
                if (param.BSIM3rdsw < 0.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: Rdsw = {0} is negative. Set to zero.\n", param.BSIM3rdsw));
                    param.BSIM3rdsw = 0.0;
                    param.BSIM3rds0 = 0.0;
                }
                if (param.BSIM3rds0 < 0.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: Rds at current temperature = {0} is negative. Set to zero.\n", param.BSIM3rds0));
                    param.BSIM3rds0 = 0.0;
                }
                if (param.BSIM3vsattemp < 1.0e3)
                    CircuitWarning.Warning(this, string.Format("Warning: Vsat at current temperature = {0} may be too small.\n", param.BSIM3vsattemp));
                if (param.BSIM3pdibl1 < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Pdibl1 = {0} is negative.\n", param.BSIM3pdibl1));
                if (param.BSIM3pdibl2 < 0.0)
                    CircuitWarning.Warning(this, string.Format("Warning: Pdibl2 = {0} is negative.\n", param.BSIM3pdibl2));
                //  Check overlap capacitance parameters 
                if (model.BSIM3cgdo < 0.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: cgdo = {0} is negative. Set to zero.\n", model.BSIM3cgdo));
                    model.BSIM3cgdo.Value = 0.0;
                }
                if (model.BSIM3cgso < 0.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: cgso = {0} is negative. Set to zero.\n", model.BSIM3cgso));
                    model.BSIM3cgso.Value = 0.0;
                }
                if (model.BSIM3cgbo < 0.0)
                {
                    CircuitWarning.Warning(this, string.Format("Warning: cgbo = {0} is negative. Set to zero.\n", model.BSIM3cgbo));
                    model.BSIM3cgbo.Value = 0.0;
                }
            }
        }

        /// <summary>
        /// Load the BSIM3 model in the circuit
        /// </summary>
        /// <param name="ckt"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            BSIM3Model model = mod as BSIM3Model;
            if (!(p is BSIM3Model.BSIM3LoadParameters))
                throw new CircuitException("No correct model parameters received. Cannot load.");
            BSIM3Model.BSIM3LoadParameters lp = (BSIM3Model.BSIM3LoadParameters)p;
            CircuitState state = ckt.State;
            var integ = ckt.Integration;

            double SourceSatCurrent, DrainSatCurrent, ag0, qgd, qgs, qgb, von, cbhat, VgstNVt, ExpVgst, cdrain, cdhat, 
                cdreq, ceqbd, ceqbs, ceqqb, ceqqd, ceqqg, czbd, czbdsw, czbdswg, czbs, czbssw, czbsswg, evbd, evbs, 
                arg, sarg, delvbd, delvbs, delvds, delvgd, delvgs, Vfbeff, dVfbeff_dVg, dVfbeff_dVb, V3, V4, gcbdb, 
                gcbgb, gcbsb, gcddb, gcdgb, gcdsb, gcgdb, gcggb, gcgsb, gcsdb, gcsgb, gcssb, MJ, MJSW, MJSWG, vbd, vbs, 
                vds, vgb, vgd, vgs, vgdo, qgate = 0.0, qbulk = 0.0, qdrn = 0.0, qsrc = 0.0, qinoi, cqgate, cqbulk, cqdrn, 
                Vds, Vgs, Vbs, Gmbs, FwdSum, RevSum, Vgs_eff, Vfb, Phis, dPhis_dVb, sqrtPhis, dsqrtPhis_dVb, Vth, 
                dVth_dVb, dVth_dVd, Vgst, dVgst_dVg, dVgst_dVb, dVgs_eff_dVg, Nvtm, Vtm, n, dn_dVb, dn_dVd, voffcv, noff, 
                dnoff_dVd, dnoff_dVb, ExpArg, V0, CoxWLcen, QovCox, LINK, DeltaPhi, dDeltaPhi_dVg, VgDP, dVgDP_dVg, Cox, 
                Tox, Tcen, dTcen_dVg, dTcen_dVd, dTcen_dVb, Ccen, Coxeff, dCoxeff_dVg, dCoxeff_dVd, dCoxeff_dVb, Denomi, 
                dDenomi_dVg, dDenomi_dVd, dDenomi_dVb, ueff, dueff_dVg, dueff_dVd, dueff_dVb, Esat, Vdsat, EsatL, 
                dEsatL_dVg, dEsatL_dVd, dEsatL_dVb, dVdsat_dVg, dVdsat_dVb, dVdsat_dVd, Vasat, dAlphaz_dVg, dAlphaz_dVb, 
                dVasat_dVg, dVasat_dVb, dVasat_dVd, Va, dVa_dVd, dVa_dVg, dVa_dVb, Vbseff, dVbseff_dVb, VbseffCV, 
                dVbseffCV_dVb, Arg1, One_Third_CoxWL, Two_Third_CoxWL, Alphaz, CoxWL, T0, dT0_dVg, dT0_dVd, dT0_dVb, 
                T1, dT1_dVg, dT1_dVd, dT1_dVb, T2, dT2_dVg, dT2_dVd, dT2_dVb, T3, dT3_dVg, dT3_dVd, dT3_dVb, T4, T5, T6, 
                T7, T8, T9, T10, T11, T12, tmp, Abulk, dAbulk_dVb, Abulk0, dAbulk0_dVb, VACLM, dVACLM_dVg, dVACLM_dVd, 
                dVACLM_dVb, VADIBL, dVADIBL_dVg, dVADIBL_dVd, dVADIBL_dVb, Xdep, dXdep_dVb, lt1, dlt1_dVb, ltw, dltw_dVb, 
                Delt_vth, dDelt_vth_dVb, Theta0, dTheta0_dVb, TempRatio, tmp1, tmp2, tmp3, tmp4, DIBL_Sft, dDIBL_Sft_dVd, 
                Lambda, dLambda_dVg, Idtot, Ibtot, a1, Vgsteff, dVgsteff_dVg, dVgsteff_dVd, dVgsteff_dVb, Vdseff, 
                dVdseff_dVg, dVdseff_dVd, dVdseff_dVb, VdseffCV, dVdseffCV_dVg, dVdseffCV_dVd, dVdseffCV_dVb, diffVds, 
                dAbulk_dVg, beta, dbeta_dVg, dbeta_dVd, dbeta_dVb, gche, dgche_dVg, dgche_dVd, dgche_dVb, fgche1, 
                dfgche1_dVg, dfgche1_dVd, dfgche1_dVb, fgche2, dfgche2_dVg, dfgche2_dVd, dfgche2_dVb, Idl, dIdl_dVg, 
                dIdl_dVd, dIdl_dVb, Idsa, dIdsa_dVg, dIdsa_dVd, dIdsa_dVb, Ids, Gm, Gds, Gmb, Isub, Gbd, Gbg, Gbb, VASCBE, 
                dVASCBE_dVg, dVASCBE_dVd, dVASCBE_dVb, CoxWovL, Rds, dRds_dVg, dRds_dVb, WVCox, WVCoxRds, Vgst2Vtm, 
                VdsatCV, dVdsatCV_dVg, dVdsatCV_dVb, Leff, Weff, dWeff_dVg, dWeff_dVb, AbulkCV, dAbulkCV_dVb, qgdo, qgso, 
                cgdo, cgso, qcheq = 0.0, qdef, gqdef = 0.0, cqdef, cqcheq, gtau_diff, gtau_drift, gcqdb = 0.0, gcqsb = 0.0, 
                gcqgb = 0.0, gcqbb = 0.0, dxpart, sxpart, ggtg, ggtd, ggts, ggtb, ddxpart_dVd, ddxpart_dVg, ddxpart_dVb, 
                ddxpart_dVs, dsxpart_dVd, dsxpart_dVg, dsxpart_dVb, dsxpart_dVs, gbspsp, gbbdp, gbbsp, gbspg, gbspb, gbspdp, 
                gbdpdp, gbdpg, gbdpb, gbdpsp, Cgg, Cgd, Cgb, Cdg, Cdd, Cds, Csg, Csd, Css, Csb, Cbg, Cbd, Cbb, Cgg1, Cgb1, 
                Cgd1, Cbg1, Cbb1, Cbd1, Qac0, Qsub0, dQac0_dVg, dQac0_dVb, dQsub0_dVg, dQsub0_dVd, dQsub0_dVb;
            bool Check, ByPass;

            Check = true;
            ByPass = false;
            if ((ckt.Mode.HasFlag(Circuit.Modes.InitSmSig)))
            {
                vbs = state.States[0][BSIM3state + BSIM3vbs];
                vgs = state.States[0][BSIM3state + BSIM3vgs];
                vds = state.States[0][BSIM3state + BSIM3vds];
                qdef = state.States[0][BSIM3state + BSIM3qdef];
            }
            else if ((ckt.Mode.HasFlag(Circuit.Modes.InitTran)))
            {
                vbs = state.States[1][BSIM3state + BSIM3vbs];
                vgs = state.States[1][BSIM3state + BSIM3vgs];
                vds = state.States[1][BSIM3state + BSIM3vds];
                qdef = state.States[1][BSIM3state + BSIM3qdef];
            }
            else if ((ckt.Mode.HasFlag(Circuit.Modes.InitJct)) && !BSIM3off)
            {
                vds = model.BSIM3type * BSIM3icVDS;
                vgs = model.BSIM3type * BSIM3icVGS;
                vbs = model.BSIM3type * BSIM3icVBS;
                qdef = 0.0;

                if ((vds == 0.0) && (vgs == 0.0) && (vbs == 0.0) &&
                        (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.DcOp |
                                          Circuit.Modes.DcTranCurve)) != 0) || (!(ckt.Mode.HasFlag(Circuit.Modes.Uic)))))
                {
                    vbs = 0.0;
                    vgs = model.BSIM3type * param.BSIM3vth0 + 0.1;
                    vds = 0.1;
                }
            }
            else if (((ckt.Mode & (Circuit.Modes.InitJct | Circuit.Modes.InitFix)) != 0) &&
                     (BSIM3off))
            {
                qdef = vbs = vgs = vds = 0.0;
            }
            else
            {
                vbs = model.BSIM3type
                      * (state.OldSolution[BSIM3bNode]
                         - state.OldSolution[BSIM3sNodePrime]);
                vgs = model.BSIM3type
                      * (state.OldSolution[BSIM3gNode]
                         - state.OldSolution[BSIM3sNodePrime]);
                vds = model.BSIM3type
                      * (state.OldSolution[BSIM3dNodePrime]
                         - state.OldSolution[BSIM3sNodePrime]);
                qdef = model.BSIM3type
                       * (state.OldSolution[BSIM3qNode]);
                vbd = vbs - vds;
                vgd = vgs - vds;
                vgdo = state.States[0][BSIM3state + BSIM3vgs]
                       - state.States[0][BSIM3state + BSIM3vds];
                delvbs = vbs - state.States[0][BSIM3state + BSIM3vbs];
                delvbd = vbd - state.States[0][BSIM3state + BSIM3vbd];
                delvgs = vgs - state.States[0][BSIM3state + BSIM3vgs];
                delvds = vds - state.States[0][BSIM3state + BSIM3vds];
                delvgd = vgd - vgdo;

                if (BSIM3mode >= 0)
                {
                    Idtot = BSIM3cd + BSIM3csub - BSIM3cbd;
                    cdhat = Idtot - BSIM3gbd * delvbd
                            + (BSIM3gmbs + BSIM3gbbs) * delvbs
                            + (BSIM3gm + BSIM3gbgs) * delvgs
                            + (BSIM3gds + BSIM3gbds) * delvds;
                    Ibtot = BSIM3cbs + BSIM3cbd - BSIM3csub;
                    cbhat = Ibtot + BSIM3gbd * delvbd
                            + (BSIM3gbs - BSIM3gbbs) * delvbs
                            - BSIM3gbgs * delvgs
                            - BSIM3gbds * delvds;
                }
                else
                {
                    Idtot = BSIM3cd - BSIM3cbd;
                    cdhat = Idtot - (BSIM3gbd - BSIM3gmbs) * delvbd
                            + BSIM3gm * delvgd
                            - BSIM3gds * delvds;
                    Ibtot = BSIM3cbs + BSIM3cbd - BSIM3csub;
                    cbhat = Ibtot + BSIM3gbs * delvbs
                            + (BSIM3gbd - BSIM3gbbs) * delvbd
                            - BSIM3gbgs * delvgd
                            + BSIM3gbds * delvds;
                }


                von = BSIM3von;
                if (state.States[0][BSIM3state + BSIM3vds] >= 0.0)
                {
                    vgs = fetlim(vgs, state.States[0][BSIM3state + BSIM3vgs], von);
                    vds = vgs - vgd;
                    vds = limvds(vds, state.States[0][BSIM3state + BSIM3vds]);
                    vgd = vgs - vds;

                }
                else
                {
                    vgd = fetlim(vgd, vgdo, von);
                    vds = vgs - vgd;
                    vds = -limvds(-vds, -(state.States[0][BSIM3state + BSIM3vds]));
                    vgs = vgd + vds;
                }

                if (vds >= 0.0)
                {
                    vbs = pnjlim(vbs, state.States[0][BSIM3state + BSIM3vbs],
                                    Circuit.CONSTvt0, model.BSIM3vcrit, ref Check);
                    vbd = vbs - vds;
                }
                else
                {
                    vbd = pnjlim(vbd, state.States[0][BSIM3state + BSIM3vbd],
                                    Circuit.CONSTvt0, model.BSIM3vcrit, ref Check);
                    vbs = vbd + vds;
                }
            }

            /* determine DC current and derivatives */
            vbd = vbs - vds;
            vgd = vgs - vds;
            vgb = vgs - vbs;

            /* Source/drain junction diode DC model begins */
            Nvtm = model.BSIM3vtm * model.BSIM3jctEmissionCoeff;
            if ((BSIM3sourceArea <= 0.0) && (BSIM3sourcePerimeter <= 0.0))
            {
                SourceSatCurrent = 1.0e-14;
            }
            else
            {
                SourceSatCurrent = BSIM3sourceArea
                                   * model.BSIM3jctTempSatCurDensity
                                   + BSIM3sourcePerimeter
                                   * model.BSIM3jctSidewallTempSatCurDensity;
            }
            if (SourceSatCurrent <= 0.0)
            {
                BSIM3gbs = ckt.Config.Gmin;
                BSIM3cbs = BSIM3gbs * vbs;
            }
            else
            {
                if (model.BSIM3ijth == 0.0)
                {
                    evbs = Math.Exp(vbs / Nvtm);
                    BSIM3gbs = SourceSatCurrent * evbs / Nvtm + ckt.Config.Gmin;
                    BSIM3cbs = SourceSatCurrent * (evbs - 1.0)
                                     + ckt.Config.Gmin * vbs;
                }
                else
                {
                    if (vbs < BSIM3vjsm)
                    {
                        evbs = Math.Exp(vbs / Nvtm);
                        BSIM3gbs = SourceSatCurrent * evbs / Nvtm + ckt.Config.Gmin;
                        BSIM3cbs = SourceSatCurrent * (evbs - 1.0)
                                         + ckt.Config.Gmin * vbs;
                    }
                    else
                    {
                        T0 = BSIM3IsEvjsm / Nvtm;
                        BSIM3gbs = T0 + ckt.Config.Gmin;
                        BSIM3cbs = BSIM3IsEvjsm - SourceSatCurrent
                                         + T0 * (vbs - BSIM3vjsm)
                                         + ckt.Config.Gmin * vbs;
                    }
                }
            }

            if ((BSIM3drainArea <= 0.0) && (BSIM3drainPerimeter <= 0.0))
            {
                DrainSatCurrent = 1.0e-14;
            }
            else
            {
                DrainSatCurrent = BSIM3drainArea
                                  * model.BSIM3jctTempSatCurDensity
                                  + BSIM3drainPerimeter
                                  * model.BSIM3jctSidewallTempSatCurDensity;
            }
            if (DrainSatCurrent <= 0.0)
            {
                BSIM3gbd = ckt.Config.Gmin;
                BSIM3cbd = BSIM3gbd * vbd;
            }
            else
            {
                if (model.BSIM3ijth == 0.0)
                {
                    evbd = Math.Exp(vbd / Nvtm);
                    BSIM3gbd = DrainSatCurrent * evbd / Nvtm + ckt.Config.Gmin;
                    BSIM3cbd = DrainSatCurrent * (evbd - 1.0)
                                     + ckt.Config.Gmin * vbd;
                }
                else
                {
                    if (vbd < BSIM3vjdm)
                    {
                        evbd = Math.Exp(vbd / Nvtm);
                        BSIM3gbd = DrainSatCurrent * evbd / Nvtm + ckt.Config.Gmin;
                        BSIM3cbd = DrainSatCurrent * (evbd - 1.0)
                                         + ckt.Config.Gmin * vbd;
                    }
                    else
                    {
                        T0 = BSIM3IsEvjdm / Nvtm;
                        BSIM3gbd = T0 + ckt.Config.Gmin;
                        BSIM3cbd = BSIM3IsEvjdm - DrainSatCurrent
                                         + T0 * (vbd - BSIM3vjdm)
                                         + ckt.Config.Gmin * vbd;
                    }
                }
            }
            /* End of diode DC model */

            if (vds >= 0.0)
            {   /* normal mode */
                BSIM3mode = 1;
                Vds = vds;
                Vgs = vgs;
                Vbs = vbs;
            }
            else
            {   /* inverse mode */
                BSIM3mode = -1;
                Vds = -vds;
                Vgs = vgd;
                Vbs = vbd;
            }

            T0 = Vbs - param.BSIM3vbsc - 0.001;
            T1 = Math.Sqrt(T0 * T0 - 0.004 * param.BSIM3vbsc);
            Vbseff = param.BSIM3vbsc + 0.5 * (T0 + T1);
            dVbseff_dVb = 0.5 * (1.0 + T0 / T1);
            if (Vbseff < Vbs)
            {
                Vbseff = Vbs;
            }

            if (Vbseff > 0.0)
            {
                T0 = param.BSIM3phi / (param.BSIM3phi + Vbseff);
                Phis = param.BSIM3phi * T0;
                dPhis_dVb = -T0 * T0;
                sqrtPhis = param.BSIM3phis3 / (param.BSIM3phi + 0.5 * Vbseff);
                dsqrtPhis_dVb = -0.5 * sqrtPhis * sqrtPhis / param.BSIM3phis3;
            }
            else
            {
                Phis = param.BSIM3phi - Vbseff;
                dPhis_dVb = -1.0;
                sqrtPhis = Math.Sqrt(Phis);
                dsqrtPhis_dVb = -0.5 / sqrtPhis;
            }
            Xdep = param.BSIM3Xdep0 * sqrtPhis / param.BSIM3sqrtPhi;
            dXdep_dVb = (param.BSIM3Xdep0 / param.BSIM3sqrtPhi)
                        * dsqrtPhis_dVb;

            Leff = param.BSIM3leff;
            Vtm = model.BSIM3vtm;
            /* Vth Calculation */
            T3 = Math.Sqrt(Xdep);
            V0 = param.BSIM3vbi - param.BSIM3phi;

            T0 = param.BSIM3dvt2 * Vbseff;
            if (T0 >= -0.5)
            {
                T1 = 1.0 + T0;
                T2 = param.BSIM3dvt2;
            }
            else /* Added to avoid any discontinuity problems caused by dvt2 */
            {
                T4 = 1.0 / (3.0 + 8.0 * T0);
                T1 = (1.0 + 3.0 * T0) * T4;
                T2 = param.BSIM3dvt2 * T4 * T4;
            }
            lt1 = model.BSIM3factor1 * T3 * T1;
            dlt1_dVb = model.BSIM3factor1 * (0.5 / T3 * T1 * dXdep_dVb + T3 * T2);

            T0 = param.BSIM3dvt2w * Vbseff;
            if (T0 >= -0.5)
            {
                T1 = 1.0 + T0;
                T2 = param.BSIM3dvt2w;
            }
            else /* Added to avoid any discontinuity problems caused by dvt2w */
            {
                T4 = 1.0 / (3.0 + 8.0 * T0);
                T1 = (1.0 + 3.0 * T0) * T4;
                T2 = param.BSIM3dvt2w * T4 * T4;
            }
            ltw = model.BSIM3factor1 * T3 * T1;
            dltw_dVb = model.BSIM3factor1 * (0.5 / T3 * T1 * dXdep_dVb + T3 * T2);

            T0 = -0.5 * param.BSIM3dvt1 * Leff / lt1;
            if (T0 > -BSIM3Model.EXP_THRESHOLD)
            {
                T1 = Math.Exp(T0);
                Theta0 = T1 * (1.0 + 2.0 * T1);
                dT1_dVb = -T0 / lt1 * T1 * dlt1_dVb;
                dTheta0_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
            }
            else
            {
                T1 = BSIM3Model.MIN_EXP;
                Theta0 = T1 * (1.0 + 2.0 * T1);
                dTheta0_dVb = 0.0;
            }

            BSIM3thetavth = param.BSIM3dvt0 * Theta0;
            Delt_vth = BSIM3thetavth * V0;
            dDelt_vth_dVb = param.BSIM3dvt0 * dTheta0_dVb * V0;

            T0 = -0.5 * param.BSIM3dvt1w * param.BSIM3weff * Leff / ltw;
            if (T0 > -BSIM3Model.EXP_THRESHOLD)
            {
                T1 = Math.Exp(T0);
                T2 = T1 * (1.0 + 2.0 * T1);
                dT1_dVb = -T0 / ltw * T1 * dltw_dVb;
                dT2_dVb = (1.0 + 4.0 * T1) * dT1_dVb;
            }
            else
            {
                T1 = BSIM3Model.MIN_EXP;
                T2 = T1 * (1.0 + 2.0 * T1);
                dT2_dVb = 0.0;
            }

            T0 = param.BSIM3dvt0w * T2;
            T2 = T0 * V0;
            dT2_dVb = param.BSIM3dvt0w * dT2_dVb * V0;

            TempRatio = ckt.Config.Temperature / model.BSIM3tnom - 1.0;
            T0 = Math.Sqrt(1.0 + param.BSIM3nlx / Leff);
            T1 = param.BSIM3k1ox * (T0 - 1.0) * param.BSIM3sqrtPhi
                 + (param.BSIM3kt1 + param.BSIM3kt1l / Leff
                    + param.BSIM3kt2 * Vbseff) * TempRatio;
            tmp2 = model.BSIM3tox * param.BSIM3phi
                   / (param.BSIM3weff + param.BSIM3w0);

            T3 = param.BSIM3eta0 + param.BSIM3etab * Vbseff;
            if (T3 < 1.0e-4) /* avoid  discontinuity problems caused by etab */
            {
                T9 = 1.0 / (3.0 - 2.0e4 * T3);
                T3 = (2.0e-4 - T3) * T9;
                T4 = T9 * T9;
            }
            else
            {
                T4 = 1.0;
            }
            dDIBL_Sft_dVd = T3 * param.BSIM3theta0vb0;
            DIBL_Sft = dDIBL_Sft_dVd * Vds;

            Vth = model.BSIM3type * param.BSIM3vth0 - param.BSIM3k1
                  * param.BSIM3sqrtPhi + param.BSIM3k1ox * sqrtPhis
                  - param.BSIM3k2ox * Vbseff - Delt_vth - T2 + (param.BSIM3k3
                          + param.BSIM3k3b * Vbseff) * tmp2 + T1 - DIBL_Sft;

            BSIM3von = Vth;

            dVth_dVb = param.BSIM3k1ox * dsqrtPhis_dVb - param.BSIM3k2ox
                       - dDelt_vth_dVb - dT2_dVb + param.BSIM3k3b * tmp2
                       - param.BSIM3etab * Vds * param.BSIM3theta0vb0 * T4
                       + param.BSIM3kt2 * TempRatio;
            dVth_dVd = -dDIBL_Sft_dVd;

            /* Calculate n */
            tmp2 = param.BSIM3nfactor * BSIM3Model.EPSSI / Xdep;
            tmp3 = param.BSIM3cdsc + param.BSIM3cdscb * Vbseff
                   + param.BSIM3cdscd * Vds;
            tmp4 = (tmp2 + tmp3 * Theta0 + param.BSIM3cit) / model.BSIM3cox;
            if (tmp4 >= -0.5)
            {
                n = 1.0 + tmp4;
                dn_dVb = (-tmp2 / Xdep * dXdep_dVb + tmp3 * dTheta0_dVb
                          + param.BSIM3cdscb * Theta0) / model.BSIM3cox;
                dn_dVd = param.BSIM3cdscd * Theta0 / model.BSIM3cox;
            }
            else
            /* avoid  discontinuity problems caused by tmp4 */
            {
                T0 = 1.0 / (3.0 + 8.0 * tmp4);
                n = (1.0 + 3.0 * tmp4) * T0;
                T0 *= T0;
                dn_dVb = (-tmp2 / Xdep * dXdep_dVb + tmp3 * dTheta0_dVb
                          + param.BSIM3cdscb * Theta0) / model.BSIM3cox * T0;
                dn_dVd = param.BSIM3cdscd * Theta0 / model.BSIM3cox * T0;
            }

            /* Poly Gate Si Depletion Effect */
            T0 = param.BSIM3vfb + param.BSIM3phi;
            if ((param.BSIM3ngate > 1.0e18) && (param.BSIM3ngate < 1.0e25)
                    && (Vgs > T0))
            /* added to avoid the problem caused by ngate */
            {
                T1 = 1.0e6 * BSIM3Model.Charge_q * BSIM3Model.EPSSI * param.BSIM3ngate
                     / (model.BSIM3cox * model.BSIM3cox);
                T4 = Math.Sqrt(1.0 + 2.0 * (Vgs - T0) / T1);
                T2 = T1 * (T4 - 1.0);
                T3 = 0.5 * T2 * T2 / T1; /* T3 = Vpoly */
                T7 = 1.12 - T3 - 0.05;
                T6 = Math.Sqrt(T7 * T7 + 0.224);
                T5 = 1.12 - 0.5 * (T7 + T6);
                Vgs_eff = Vgs - T5;
                dVgs_eff_dVg = 1.0 - (0.5 - 0.5 / T4) * (1.0 + T7 / T6);
            }
            else
            {
                Vgs_eff = Vgs;
                dVgs_eff_dVg = 1.0;
            }
            Vgst = Vgs_eff - Vth;

            /* Effective Vgst (Vgsteff) Calculation */

            T10 = 2.0 * n * Vtm;
            VgstNVt = Vgst / T10;
            ExpArg = (2.0 * param.BSIM3voff - Vgst) / T10;

            /* MCJ: Very small Vgst */
            if (VgstNVt > BSIM3Model.EXP_THRESHOLD)
            {
                Vgsteff = Vgst;
                dVgsteff_dVg = dVgs_eff_dVg;
                dVgsteff_dVd = -dVth_dVd;
                dVgsteff_dVb = -dVth_dVb;
            }
            else if (ExpArg > BSIM3Model.EXP_THRESHOLD)
            {
                T0 = (Vgst - param.BSIM3voff) / (n * Vtm);
                ExpVgst = Math.Exp(T0);
                Vgsteff = Vtm * param.BSIM3cdep0 / model.BSIM3cox * ExpVgst;
                dVgsteff_dVg = Vgsteff / (n * Vtm);
                dVgsteff_dVd = -dVgsteff_dVg * (dVth_dVd + T0 * Vtm * dn_dVd);
                dVgsteff_dVb = -dVgsteff_dVg * (dVth_dVb + T0 * Vtm * dn_dVb);
                dVgsteff_dVg *= dVgs_eff_dVg;
            }
            else
            {
                ExpVgst = Math.Exp(VgstNVt);
                T1 = T10 * Math.Log(1.0 + ExpVgst);
                dT1_dVg = ExpVgst / (1.0 + ExpVgst);
                dT1_dVb = -dT1_dVg * (dVth_dVb + Vgst / n * dn_dVb)
                          + T1 / n * dn_dVb;
                dT1_dVd = -dT1_dVg * (dVth_dVd + Vgst / n * dn_dVd)
                          + T1 / n * dn_dVd;

                dT2_dVg = -model.BSIM3cox / (Vtm * param.BSIM3cdep0)
                          * Math.Exp(ExpArg);
                T2 = 1.0 - T10 * dT2_dVg;
                dT2_dVd = -dT2_dVg * (dVth_dVd - 2.0 * Vtm * ExpArg * dn_dVd)
                          + (T2 - 1.0) / n * dn_dVd;
                dT2_dVb = -dT2_dVg * (dVth_dVb - 2.0 * Vtm * ExpArg * dn_dVb)
                          + (T2 - 1.0) / n * dn_dVb;

                Vgsteff = T1 / T2;
                T3 = T2 * T2;
                dVgsteff_dVg = (T2 * dT1_dVg - T1 * dT2_dVg) / T3 * dVgs_eff_dVg;
                dVgsteff_dVd = (T2 * dT1_dVd - T1 * dT2_dVd) / T3;
                dVgsteff_dVb = (T2 * dT1_dVb - T1 * dT2_dVb) / T3;
            }
            BSIM3Vgsteff = Vgsteff;

            /* Calculate Effective Channel Geometry */
            T9 = sqrtPhis - param.BSIM3sqrtPhi;
            Weff = param.BSIM3weff - 2.0 * (param.BSIM3dwg * Vgsteff
                                              + param.BSIM3dwb * T9);
            dWeff_dVg = -2.0 * param.BSIM3dwg;
            dWeff_dVb = -2.0 * param.BSIM3dwb * dsqrtPhis_dVb;

            if (Weff < 2.0e-8) /* to avoid the discontinuity problem due to Weff*/
            {
                T0 = 1.0 / (6.0e-8 - 2.0 * Weff);
                Weff = 2.0e-8 * (4.0e-8 - Weff) * T0;
                T0 *= T0 * 4.0e-16;
                dWeff_dVg *= T0;
                dWeff_dVb *= T0;
            }

            T0 = param.BSIM3prwg * Vgsteff + param.BSIM3prwb * T9;
            if (T0 >= -0.9)
            {
                Rds = param.BSIM3rds0 * (1.0 + T0);
                dRds_dVg = param.BSIM3rds0 * param.BSIM3prwg;
                dRds_dVb = param.BSIM3rds0 * param.BSIM3prwb * dsqrtPhis_dVb;
            }
            else
            /* to avoid the discontinuity problem due to prwg and prwb*/
            {
                T1 = 1.0 / (17.0 + 20.0 * T0);
                Rds = param.BSIM3rds0 * (0.8 + T0) * T1;
                T1 *= T1;
                dRds_dVg = param.BSIM3rds0 * param.BSIM3prwg * T1;
                dRds_dVb = param.BSIM3rds0 * param.BSIM3prwb * dsqrtPhis_dVb
                           * T1;
            }
            BSIM3rds = Rds; /* Noise Bugfix */

            /* Calculate Abulk */
            T1 = 0.5 * param.BSIM3k1ox / sqrtPhis;
            dT1_dVb = -T1 / sqrtPhis * dsqrtPhis_dVb;

            T9 = Math.Sqrt(param.BSIM3xj * Xdep);
            tmp1 = Leff + 2.0 * T9;
            T5 = Leff / tmp1;
            tmp2 = param.BSIM3a0 * T5;
            tmp3 = param.BSIM3weff + param.BSIM3b1;
            tmp4 = param.BSIM3b0 / tmp3;
            T2 = tmp2 + tmp4;
            dT2_dVb = -T9 / tmp1 / Xdep * dXdep_dVb;
            T6 = T5 * T5;
            T7 = T5 * T6;

            Abulk0 = 1.0 + T1 * T2;
            dAbulk0_dVb = T1 * tmp2 * dT2_dVb + T2 * dT1_dVb;

            T8 = param.BSIM3ags * param.BSIM3a0 * T7;
            dAbulk_dVg = -T1 * T8;
            Abulk = Abulk0 + dAbulk_dVg * Vgsteff;
            dAbulk_dVb = dAbulk0_dVb - T8 * Vgsteff * (dT1_dVb
                         + 3.0 * T1 * dT2_dVb);

            if (Abulk0 < 0.1) /* added to avoid the problems caused by Abulk0 */
            {
                T9 = 1.0 / (3.0 - 20.0 * Abulk0);
                Abulk0 = (0.2 - Abulk0) * T9;
                dAbulk0_dVb *= T9 * T9;
            }

            if (Abulk < 0.1)
            /* added to avoid the problems caused by Abulk */
            {
                T9 = 1.0 / (3.0 - 20.0 * Abulk);
                Abulk = (0.2 - Abulk) * T9;
                T10 = T9 * T9;
                dAbulk_dVb *= T10;
                dAbulk_dVg *= T10;
            }
            BSIM3Abulk = Abulk;

            T2 = param.BSIM3keta * Vbseff;
            if (T2 >= -0.9)
            {
                T0 = 1.0 / (1.0 + T2);
                dT0_dVb = -param.BSIM3keta * T0 * T0;
            }
            else
            /* added to avoid the problems caused by Keta */
            {
                T1 = 1.0 / (0.8 + T2);
                T0 = (17.0 + 20.0 * T2) * T1;
                dT0_dVb = -param.BSIM3keta * T1 * T1;
            }
            dAbulk_dVg *= T0;
            dAbulk_dVb = dAbulk_dVb * T0 + Abulk * dT0_dVb;
            dAbulk0_dVb = dAbulk0_dVb * T0 + Abulk0 * dT0_dVb;
            Abulk *= T0;
            Abulk0 *= T0;


            /* Mobility calculation */
            if (model.BSIM3mobMod == 1)
            {
                T0 = Vgsteff + Vth + Vth;
                T2 = param.BSIM3ua + param.BSIM3uc * Vbseff;
                T3 = T0 / model.BSIM3tox;
                T5 = T3 * (T2 + param.BSIM3ub * T3);
                dDenomi_dVg = (T2 + 2.0 * param.BSIM3ub * T3) / model.BSIM3tox;
                dDenomi_dVd = dDenomi_dVg * 2.0 * dVth_dVd;
                dDenomi_dVb = dDenomi_dVg * 2.0 * dVth_dVb + param.BSIM3uc * T3;
            }
            else if (model.BSIM3mobMod == 2)
            {
                T5 = Vgsteff / model.BSIM3tox * (param.BSIM3ua
                                                  + param.BSIM3uc * Vbseff + param.BSIM3ub * Vgsteff
                                                  / model.BSIM3tox);
                dDenomi_dVg = (param.BSIM3ua + param.BSIM3uc * Vbseff
                               + 2.0 * param.BSIM3ub * Vgsteff / model.BSIM3tox)
                              / model.BSIM3tox;
                dDenomi_dVd = 0.0;
                dDenomi_dVb = Vgsteff * param.BSIM3uc / model.BSIM3tox;
            }
            else
            {
                T0 = Vgsteff + Vth + Vth;
                T2 = 1.0 + param.BSIM3uc * Vbseff;
                T3 = T0 / model.BSIM3tox;
                T4 = T3 * (param.BSIM3ua + param.BSIM3ub * T3);
                T5 = T4 * T2;
                dDenomi_dVg = (param.BSIM3ua + 2.0 * param.BSIM3ub * T3) * T2
                              / model.BSIM3tox;
                dDenomi_dVd = dDenomi_dVg * 2.0 * dVth_dVd;
                dDenomi_dVb = dDenomi_dVg * 2.0 * dVth_dVb + param.BSIM3uc * T4;
            }

            if (T5 >= -0.8)
            {
                Denomi = 1.0 + T5;
            }
            else /* Added to avoid the discontinuity problem caused by ua and ub*/
            {
                T9 = 1.0 / (7.0 + 10.0 * T5);
                Denomi = (0.6 + T5) * T9;
                T9 *= T9;
                dDenomi_dVg *= T9;
                dDenomi_dVd *= T9;
                dDenomi_dVb *= T9;
            }

            BSIM3ueff = ueff = param.BSIM3u0temp / Denomi;
            T9 = -ueff / Denomi;
            dueff_dVg = T9 * dDenomi_dVg;
            dueff_dVd = T9 * dDenomi_dVd;
            dueff_dVb = T9 * dDenomi_dVb;

            /* Saturation Drain Voltage  Vdsat */
            WVCox = Weff * param.BSIM3vsattemp * model.BSIM3cox;
            WVCoxRds = WVCox * Rds;

            Esat = 2.0 * param.BSIM3vsattemp / ueff;
            EsatL = Esat * Leff;
            T0 = -EsatL / ueff;
            dEsatL_dVg = T0 * dueff_dVg;
            dEsatL_dVd = T0 * dueff_dVd;
            dEsatL_dVb = T0 * dueff_dVb;

            /* Sqrt() */
            a1 = param.BSIM3a1;
            if (a1 == 0.0)
            {
                Lambda = param.BSIM3a2;
                dLambda_dVg = 0.0;
            }
            else if (a1 > 0.0)
            /* Added to avoid the discontinuity problem
               caused by a1 and a2 (Lambda) */
            {
                T0 = 1.0 - param.BSIM3a2;
                T1 = T0 - param.BSIM3a1 * Vgsteff - 0.0001;
                T2 = Math.Sqrt(T1 * T1 + 0.0004 * T0);
                Lambda = param.BSIM3a2 + T0 - 0.5 * (T1 + T2);
                dLambda_dVg = 0.5 * param.BSIM3a1 * (1.0 + T1 / T2);
            }
            else
            {
                T1 = param.BSIM3a2 + param.BSIM3a1 * Vgsteff - 0.0001;
                T2 = Math.Sqrt(T1 * T1 + 0.0004 * param.BSIM3a2);
                Lambda = 0.5 * (T1 + T2);
                dLambda_dVg = 0.5 * param.BSIM3a1 * (1.0 + T1 / T2);
            }

            Vgst2Vtm = Vgsteff + 2.0 * Vtm;
            BSIM3AbovVgst2Vtm = Abulk / Vgst2Vtm;

            if (Rds > 0)
            {
                tmp2 = dRds_dVg / Rds + dWeff_dVg / Weff;
                tmp3 = dRds_dVb / Rds + dWeff_dVb / Weff;
            }
            else
            {
                tmp2 = dWeff_dVg / Weff;
                tmp3 = dWeff_dVb / Weff;
            }
            if ((Rds == 0.0) && (Lambda == 1.0))
            {
                T0 = 1.0 / (Abulk * EsatL + Vgst2Vtm);
                tmp1 = 0.0;
                T1 = T0 * T0;
                T2 = Vgst2Vtm * T0;
                T3 = EsatL * Vgst2Vtm;
                Vdsat = T3 * T0;

                dT0_dVg = -(Abulk * dEsatL_dVg + EsatL * dAbulk_dVg + 1.0) * T1;
                dT0_dVd = -(Abulk * dEsatL_dVd) * T1;
                dT0_dVb = -(Abulk * dEsatL_dVb + dAbulk_dVb * EsatL) * T1;

                dVdsat_dVg = T3 * dT0_dVg + T2 * dEsatL_dVg + EsatL * T0;
                dVdsat_dVd = T3 * dT0_dVd + T2 * dEsatL_dVd;
                dVdsat_dVb = T3 * dT0_dVb + T2 * dEsatL_dVb;
            }
            else
            {
                tmp1 = dLambda_dVg / (Lambda * Lambda);
                T9 = Abulk * WVCoxRds;
                T8 = Abulk * T9;
                T7 = Vgst2Vtm * T9;
                T6 = Vgst2Vtm * WVCoxRds;
                T0 = 2.0 * Abulk * (T9 - 1.0 + 1.0 / Lambda);
                dT0_dVg = 2.0 * (T8 * tmp2 - Abulk * tmp1
                                 + (2.0 * T9 + 1.0 / Lambda - 1.0) * dAbulk_dVg);

                dT0_dVb = 2.0 * (T8 * (2.0 / Abulk * dAbulk_dVb + tmp3)
                                 + (1.0 / Lambda - 1.0) * dAbulk_dVb);
                dT0_dVd = 0.0;
                T1 = Vgst2Vtm * (2.0 / Lambda - 1.0) + Abulk * EsatL + 3.0 * T7;

                dT1_dVg = (2.0 / Lambda - 1.0) - 2.0 * Vgst2Vtm * tmp1
                          + Abulk * dEsatL_dVg + EsatL * dAbulk_dVg + 3.0 * (T9
                                  + T7 * tmp2 + T6 * dAbulk_dVg);
                dT1_dVb = Abulk * dEsatL_dVb + EsatL * dAbulk_dVb
                          + 3.0 * (T6 * dAbulk_dVb + T7 * tmp3);
                dT1_dVd = Abulk * dEsatL_dVd;

                T2 = Vgst2Vtm * (EsatL + 2.0 * T6);
                dT2_dVg = EsatL + Vgst2Vtm * dEsatL_dVg
                          + T6 * (4.0 + 2.0 * Vgst2Vtm * tmp2);
                dT2_dVb = Vgst2Vtm * (dEsatL_dVb + 2.0 * T6 * tmp3);
                dT2_dVd = Vgst2Vtm * dEsatL_dVd;

                T3 = Math.Sqrt(T1 * T1 - 2.0 * T0 * T2);
                Vdsat = (T1 - T3) / T0;

                dT3_dVg = (T1 * dT1_dVg - 2.0 * (T0 * dT2_dVg + T2 * dT0_dVg))
                          / T3;
                dT3_dVd = (T1 * dT1_dVd - 2.0 * (T0 * dT2_dVd + T2 * dT0_dVd))
                          / T3;
                dT3_dVb = (T1 * dT1_dVb - 2.0 * (T0 * dT2_dVb + T2 * dT0_dVb))
                          / T3;

                dVdsat_dVg = (dT1_dVg - (T1 * dT1_dVg - dT0_dVg * T2
                                         - T0 * dT2_dVg) / T3 - Vdsat * dT0_dVg) / T0;
                dVdsat_dVb = (dT1_dVb - (T1 * dT1_dVb - dT0_dVb * T2
                                         - T0 * dT2_dVb) / T3 - Vdsat * dT0_dVb) / T0;
                dVdsat_dVd = (dT1_dVd - (T1 * dT1_dVd - T0 * dT2_dVd) / T3) / T0;
            }
            BSIM3vdsat = Vdsat;

            /* Effective Vds (Vdseff) Calculation */
            T1 = Vdsat - Vds - param.BSIM3delta;
            dT1_dVg = dVdsat_dVg;
            dT1_dVd = dVdsat_dVd - 1.0;
            dT1_dVb = dVdsat_dVb;

            T2 = Math.Sqrt(T1 * T1 + 4.0 * param.BSIM3delta * Vdsat);
            T0 = T1 / T2;
            T3 = 2.0 * param.BSIM3delta / T2;
            dT2_dVg = T0 * dT1_dVg + T3 * dVdsat_dVg;
            dT2_dVd = T0 * dT1_dVd + T3 * dVdsat_dVd;
            dT2_dVb = T0 * dT1_dVb + T3 * dVdsat_dVb;

            Vdseff = Vdsat - 0.5 * (T1 + T2);
            dVdseff_dVg = dVdsat_dVg - 0.5 * (dT1_dVg + dT2_dVg);
            dVdseff_dVd = dVdsat_dVd - 0.5 * (dT1_dVd + dT2_dVd);
            dVdseff_dVb = dVdsat_dVb - 0.5 * (dT1_dVb + dT2_dVb);
            /* Added to eliminate non-zero Vdseff at Vds=0.0 */
            if (Vds == 0.0)
            {
                Vdseff = 0.0;
                dVdseff_dVg = 0.0;
                dVdseff_dVb = 0.0;
            }

            /* Calculate VAsat */
            tmp4 = 1.0 - 0.5 * Abulk * Vdsat / Vgst2Vtm;
            T9 = WVCoxRds * Vgsteff;
            T8 = T9 / Vgst2Vtm;
            T0 = EsatL + Vdsat + 2.0 * T9 * tmp4;

            T7 = 2.0 * WVCoxRds * tmp4;
            dT0_dVg = dEsatL_dVg + dVdsat_dVg + T7 * (1.0 + tmp2 * Vgsteff)
                      - T8 * (Abulk * dVdsat_dVg - Abulk * Vdsat / Vgst2Vtm
                              + Vdsat * dAbulk_dVg);

            dT0_dVb = dEsatL_dVb + dVdsat_dVb + T7 * tmp3 * Vgsteff
                      - T8 * (dAbulk_dVb * Vdsat + Abulk * dVdsat_dVb);
            dT0_dVd = dEsatL_dVd + dVdsat_dVd - T8 * Abulk * dVdsat_dVd;

            T9 = WVCoxRds * Abulk;
            T1 = 2.0 / Lambda - 1.0 + T9;
            dT1_dVg = -2.0 * tmp1 + WVCoxRds * (Abulk * tmp2 + dAbulk_dVg);
            dT1_dVb = dAbulk_dVb * WVCoxRds + T9 * tmp3;

            Vasat = T0 / T1;
            dVasat_dVg = (dT0_dVg - Vasat * dT1_dVg) / T1;
            dVasat_dVb = (dT0_dVb - Vasat * dT1_dVb) / T1;
            dVasat_dVd = dT0_dVd / T1;

            if (Vdseff > Vds)
                Vdseff = Vds;
            diffVds = Vds - Vdseff;
            BSIM3Vdseff = Vdseff;

            /* Calculate VACLM */
            if ((param.BSIM3pclm > 0.0) && (diffVds > 1.0e-10))
            {
                T0 = 1.0 / (param.BSIM3pclm * Abulk * param.BSIM3litl);
                dT0_dVb = -T0 / Abulk * dAbulk_dVb;
                dT0_dVg = -T0 / Abulk * dAbulk_dVg;

                T2 = Vgsteff / EsatL;
                T1 = Leff * (Abulk + T2);
                dT1_dVg = Leff * ((1.0 - T2 * dEsatL_dVg) / EsatL + dAbulk_dVg);
                dT1_dVb = Leff * (dAbulk_dVb - T2 * dEsatL_dVb / EsatL);
                dT1_dVd = -T2 * dEsatL_dVd / Esat;

                T9 = T0 * T1;
                VACLM = T9 * diffVds;
                dVACLM_dVg = T0 * dT1_dVg * diffVds - T9 * dVdseff_dVg
                             + T1 * diffVds * dT0_dVg;
                dVACLM_dVb = (dT0_dVb * T1 + T0 * dT1_dVb) * diffVds
                             - T9 * dVdseff_dVb;
                dVACLM_dVd = T0 * dT1_dVd * diffVds + T9 * (1.0 - dVdseff_dVd);
            }
            else
            {
                VACLM = BSIM3Model.MAX_EXP;
                dVACLM_dVd = dVACLM_dVg = dVACLM_dVb = 0.0;
            }

            /* Calculate VADIBL */
            if (param.BSIM3thetaRout > 0.0)
            {
                T8 = Abulk * Vdsat;
                T0 = Vgst2Vtm * T8;
                dT0_dVg = Vgst2Vtm * Abulk * dVdsat_dVg + T8
                          + Vgst2Vtm * Vdsat * dAbulk_dVg;
                dT0_dVb = Vgst2Vtm * (dAbulk_dVb * Vdsat + Abulk * dVdsat_dVb);
                dT0_dVd = Vgst2Vtm * Abulk * dVdsat_dVd;

                T1 = Vgst2Vtm + T8;
                dT1_dVg = 1.0 + Abulk * dVdsat_dVg + Vdsat * dAbulk_dVg;
                dT1_dVb = Abulk * dVdsat_dVb + dAbulk_dVb * Vdsat;
                dT1_dVd = Abulk * dVdsat_dVd;

                T9 = T1 * T1;
                T2 = param.BSIM3thetaRout;
                VADIBL = (Vgst2Vtm - T0 / T1) / T2;
                dVADIBL_dVg = (1.0 - dT0_dVg / T1 + T0 * dT1_dVg / T9) / T2;
                dVADIBL_dVb = (-dT0_dVb / T1 + T0 * dT1_dVb / T9) / T2;
                dVADIBL_dVd = (-dT0_dVd / T1 + T0 * dT1_dVd / T9) / T2;

                T7 = param.BSIM3pdiblb * Vbseff;
                if (T7 >= -0.9)
                {
                    T3 = 1.0 / (1.0 + T7);
                    VADIBL *= T3;
                    dVADIBL_dVg *= T3;
                    dVADIBL_dVb = (dVADIBL_dVb - VADIBL * param.BSIM3pdiblb)
                                  * T3;
                    dVADIBL_dVd *= T3;
                }
                else
                /* Added to avoid the discontinuity problem caused by pdiblcb */
                {
                    T4 = 1.0 / (0.8 + T7);
                    T3 = (17.0 + 20.0 * T7) * T4;
                    dVADIBL_dVg *= T3;
                    dVADIBL_dVb = dVADIBL_dVb * T3
                                  - VADIBL * param.BSIM3pdiblb * T4 * T4;
                    dVADIBL_dVd *= T3;
                    VADIBL *= T3;
                }
            }
            else
            {
                VADIBL = BSIM3Model.MAX_EXP;
                dVADIBL_dVd = dVADIBL_dVg = dVADIBL_dVb = 0.0;
            }

            /* Calculate VA */

            T8 = param.BSIM3pvag / EsatL;
            T9 = T8 * Vgsteff;
            if (T9 > -0.9)
            {
                T0 = 1.0 + T9;
                dT0_dVg = T8 * (1.0 - Vgsteff * dEsatL_dVg / EsatL);
                dT0_dVb = -T9 * dEsatL_dVb / EsatL;
                dT0_dVd = -T9 * dEsatL_dVd / EsatL;
            }
            else /* Added to avoid the discontinuity problems caused by pvag */
            {
                T1 = 1.0 / (17.0 + 20.0 * T9);
                T0 = (0.8 + T9) * T1;
                T1 *= T1;
                dT0_dVg = T8 * (1.0 - Vgsteff * dEsatL_dVg / EsatL) * T1;

                T9 *= T1 / EsatL;
                dT0_dVb = -T9 * dEsatL_dVb;
                dT0_dVd = -T9 * dEsatL_dVd;
            }

            tmp1 = VACLM * VACLM;
            tmp2 = VADIBL * VADIBL;
            tmp3 = VACLM + VADIBL;

            T1 = VACLM * VADIBL / tmp3;
            tmp3 *= tmp3;
            dT1_dVg = (tmp1 * dVADIBL_dVg + tmp2 * dVACLM_dVg) / tmp3;
            dT1_dVd = (tmp1 * dVADIBL_dVd + tmp2 * dVACLM_dVd) / tmp3;
            dT1_dVb = (tmp1 * dVADIBL_dVb + tmp2 * dVACLM_dVb) / tmp3;

            Va = Vasat + T0 * T1;
            dVa_dVg = dVasat_dVg + T1 * dT0_dVg + T0 * dT1_dVg;
            dVa_dVd = dVasat_dVd + T1 * dT0_dVd + T0 * dT1_dVd;
            dVa_dVb = dVasat_dVb + T1 * dT0_dVb + T0 * dT1_dVb;

            /* Calculate VASCBE */
            if (param.BSIM3pscbe2 > 0.0)
            {
                if (diffVds > param.BSIM3pscbe1 * param.BSIM3litl
                        / BSIM3Model.EXP_THRESHOLD)
                {
                    T0 = param.BSIM3pscbe1 * param.BSIM3litl / diffVds;
                    VASCBE = Leff * Math.Exp(T0) / param.BSIM3pscbe2;
                    T1 = T0 * VASCBE / diffVds;
                    dVASCBE_dVg = T1 * dVdseff_dVg;
                    dVASCBE_dVd = -T1 * (1.0 - dVdseff_dVd);
                    dVASCBE_dVb = T1 * dVdseff_dVb;
                }
                else
                {
                    VASCBE = BSIM3Model.MAX_EXP * Leff / param.BSIM3pscbe2;
                    dVASCBE_dVg = dVASCBE_dVd = dVASCBE_dVb = 0.0;
                }
            }
            else
            {
                VASCBE = BSIM3Model.MAX_EXP;
                dVASCBE_dVg = dVASCBE_dVd = dVASCBE_dVb = 0.0;
            }

            /* Calculate Ids */
            CoxWovL = model.BSIM3cox * Weff / Leff;
            beta = ueff * CoxWovL;
            dbeta_dVg = CoxWovL * dueff_dVg + beta * dWeff_dVg / Weff;
            dbeta_dVd = CoxWovL * dueff_dVd;
            dbeta_dVb = CoxWovL * dueff_dVb + beta * dWeff_dVb / Weff;

            T0 = 1.0 - 0.5 * Abulk * Vdseff / Vgst2Vtm;
            dT0_dVg = -0.5 * (Abulk * dVdseff_dVg
                              - Abulk * Vdseff / Vgst2Vtm + Vdseff * dAbulk_dVg) / Vgst2Vtm;
            dT0_dVd = -0.5 * Abulk * dVdseff_dVd / Vgst2Vtm;
            dT0_dVb = -0.5 * (Abulk * dVdseff_dVb + dAbulk_dVb * Vdseff)
                      / Vgst2Vtm;

            fgche1 = Vgsteff * T0;
            dfgche1_dVg = Vgsteff * dT0_dVg + T0;
            dfgche1_dVd = Vgsteff * dT0_dVd;
            dfgche1_dVb = Vgsteff * dT0_dVb;

            T9 = Vdseff / EsatL;
            fgche2 = 1.0 + T9;
            dfgche2_dVg = (dVdseff_dVg - T9 * dEsatL_dVg) / EsatL;
            dfgche2_dVd = (dVdseff_dVd - T9 * dEsatL_dVd) / EsatL;
            dfgche2_dVb = (dVdseff_dVb - T9 * dEsatL_dVb) / EsatL;

            gche = beta * fgche1 / fgche2;
            dgche_dVg = (beta * dfgche1_dVg + fgche1 * dbeta_dVg
                         - gche * dfgche2_dVg) / fgche2;
            dgche_dVd = (beta * dfgche1_dVd + fgche1 * dbeta_dVd
                         - gche * dfgche2_dVd) / fgche2;
            dgche_dVb = (beta * dfgche1_dVb + fgche1 * dbeta_dVb
                         - gche * dfgche2_dVb) / fgche2;

            T0 = 1.0 + gche * Rds;
            T9 = Vdseff / T0;
            Idl = gche * T9;

            dIdl_dVg = (gche * dVdseff_dVg + T9 * dgche_dVg) / T0
                       - Idl * gche / T0 * dRds_dVg;

            dIdl_dVd = (gche * dVdseff_dVd + T9 * dgche_dVd) / T0;
            dIdl_dVb = (gche * dVdseff_dVb + T9 * dgche_dVb
                        - Idl * dRds_dVb * gche) / T0;

            T9 = diffVds / Va;
            T0 = 1.0 + T9;
            Idsa = Idl * T0;
            dIdsa_dVg = T0 * dIdl_dVg - Idl * (dVdseff_dVg + T9 * dVa_dVg) / Va;
            dIdsa_dVd = T0 * dIdl_dVd + Idl * (1.0 - dVdseff_dVd
                                               - T9 * dVa_dVd) / Va;
            dIdsa_dVb = T0 * dIdl_dVb - Idl * (dVdseff_dVb + T9 * dVa_dVb) / Va;

            T9 = diffVds / VASCBE;
            T0 = 1.0 + T9;
            Ids = Idsa * T0;

            Gm = T0 * dIdsa_dVg - Idsa * (dVdseff_dVg + T9 * dVASCBE_dVg) / VASCBE;
            Gds = T0 * dIdsa_dVd + Idsa * (1.0 - dVdseff_dVd
                                           - T9 * dVASCBE_dVd) / VASCBE;
            Gmb = T0 * dIdsa_dVb - Idsa * (dVdseff_dVb
                                           + T9 * dVASCBE_dVb) / VASCBE;

            Gds += Gm * dVgsteff_dVd;
            Gmb += Gm * dVgsteff_dVb;
            Gm *= dVgsteff_dVg;
            Gmb *= dVbseff_dVb;

            /* Substrate current begins */
            tmp = param.BSIM3alpha0 + param.BSIM3alpha1 * Leff;
            if ((tmp <= 0.0) || (param.BSIM3beta0 <= 0.0))
            {
                Isub = Gbd = Gbb = Gbg = 0.0;
            }
            else
            {
                T2 = tmp / Leff;
                if (diffVds > param.BSIM3beta0 / BSIM3Model.EXP_THRESHOLD)
                {
                    T0 = -param.BSIM3beta0 / diffVds;
                    T1 = T2 * diffVds * Math.Exp(T0);
                    T3 = T1 / diffVds * (T0 - 1.0);
                    dT1_dVg = T3 * dVdseff_dVg;
                    dT1_dVd = T3 * (dVdseff_dVd - 1.0);
                    dT1_dVb = T3 * dVdseff_dVb;
                }
                else
                {
                    T3 = T2 * BSIM3Model.MIN_EXP;
                    T1 = T3 * diffVds;
                    dT1_dVg = -T3 * dVdseff_dVg;
                    dT1_dVd = T3 * (1.0 - dVdseff_dVd);
                    dT1_dVb = -T3 * dVdseff_dVb;
                }
                Isub = T1 * Idsa;
                Gbg = T1 * dIdsa_dVg + Idsa * dT1_dVg;
                Gbd = T1 * dIdsa_dVd + Idsa * dT1_dVd;
                Gbb = T1 * dIdsa_dVb + Idsa * dT1_dVb;

                Gbd += Gbg * dVgsteff_dVd;
                Gbb += Gbg * dVgsteff_dVb;
                Gbg *= dVgsteff_dVg;
                Gbb *= dVbseff_dVb; /* bug fixing */
            }

            cdrain = Ids;
            BSIM3gds = Gds;
            BSIM3gm = Gm;
            BSIM3gmbs = Gmb;

            BSIM3gbbs = Gbb;
            BSIM3gbgs = Gbg;
            BSIM3gbds = Gbd;

            BSIM3csub = Isub;

            /* BSIM3 thermal noise Qinv calculated from all capMod
             * 0, 1, 2 & 3 stored in BSIM3qinv 1/1998 */

            if ((model.BSIM3xpart < 0) || (!lp.ChargeComputationNeeded))
            {
                qgate = qdrn = qsrc = qbulk = 0.0;
                BSIM3cggb = BSIM3cgsb = BSIM3cgdb = 0.0;
                BSIM3cdgb = BSIM3cdsb = BSIM3cddb = 0.0;
                BSIM3cbgb = BSIM3cbsb = BSIM3cbdb = 0.0;
                BSIM3cqdb = BSIM3cqsb = BSIM3cqgb
                                                    = BSIM3cqbb = 0.0;
                BSIM3gtau = 0.0;
                goto finished;
            }
            else if (model.BSIM3capMod == 0)
            {
                if (Vbseff < 0.0)
                {
                    Vbseff = Vbs;
                    dVbseff_dVb = 1.0;
                }
                else
                {
                    Vbseff = param.BSIM3phi - Phis;
                    dVbseff_dVb = -dPhis_dVb;
                }

                Vfb = param.BSIM3vfbcv;
                Vth = Vfb + param.BSIM3phi + param.BSIM3k1ox * sqrtPhis;
                Vgst = Vgs_eff - Vth;
                dVth_dVb = param.BSIM3k1ox * dsqrtPhis_dVb;
                dVgst_dVb = -dVth_dVb;
                dVgst_dVg = dVgs_eff_dVg;

                CoxWL = model.BSIM3cox * param.BSIM3weffCV
                        * param.BSIM3leffCV;
                Arg1 = Vgs_eff - Vbseff - Vfb;

                if (Arg1 <= 0.0)
                {
                    qgate = CoxWL * Arg1;
                    qbulk = -qgate;
                    qdrn = 0.0;

                    BSIM3cggb = CoxWL * dVgs_eff_dVg;
                    BSIM3cgdb = 0.0;
                    BSIM3cgsb = CoxWL * (dVbseff_dVb - dVgs_eff_dVg);

                    BSIM3cdgb = 0.0;
                    BSIM3cddb = 0.0;
                    BSIM3cdsb = 0.0;

                    BSIM3cbgb = -CoxWL * dVgs_eff_dVg;
                    BSIM3cbdb = 0.0;
                    BSIM3cbsb = -BSIM3cgsb;
                    BSIM3qinv = 0.0;
                }
                else if (Vgst <= 0.0)
                {
                    T1 = 0.5 * param.BSIM3k1ox;
                    T2 = Math.Sqrt(T1 * T1 + Arg1);
                    qgate = CoxWL * param.BSIM3k1ox * (T2 - T1);
                    qbulk = -qgate;
                    qdrn = 0.0;

                    T0 = CoxWL * T1 / T2;
                    BSIM3cggb = T0 * dVgs_eff_dVg;
                    BSIM3cgdb = 0.0;
                    BSIM3cgsb = T0 * (dVbseff_dVb - dVgs_eff_dVg);

                    BSIM3cdgb = 0.0;
                    BSIM3cddb = 0.0;
                    BSIM3cdsb = 0.0;

                    BSIM3cbgb = -BSIM3cggb;
                    BSIM3cbdb = 0.0;
                    BSIM3cbsb = -BSIM3cgsb;
                    BSIM3qinv = 0.0;
                }
                else
                {
                    One_Third_CoxWL = CoxWL / 3.0;
                    Two_Third_CoxWL = 2.0 * One_Third_CoxWL;

                    AbulkCV = Abulk0 * param.BSIM3abulkCVfactor;
                    dAbulkCV_dVb = param.BSIM3abulkCVfactor * dAbulk0_dVb;
                    Vdsat = Vgst / AbulkCV;
                    dVdsat_dVg = dVgs_eff_dVg / AbulkCV;
                    dVdsat_dVb = -(Vdsat * dAbulkCV_dVb + dVth_dVb) / AbulkCV;

                    if (model.BSIM3xpart > 0.5)
                    {   /* 0/100 Charge partition model */
                        if (Vdsat <= Vds)
                        {   /* saturation region */
                            T1 = Vdsat / 3.0;
                            qgate = CoxWL * (Vgs_eff - Vfb
                                             - param.BSIM3phi - T1);
                            T2 = -Two_Third_CoxWL * Vgst;
                            qbulk = -(qgate + T2);
                            qdrn = 0.0;

                            BSIM3cggb = One_Third_CoxWL * (3.0
                                                                 - dVdsat_dVg) * dVgs_eff_dVg;
                            T2 = -One_Third_CoxWL * dVdsat_dVb;
                            BSIM3cgsb = -(BSIM3cggb + T2);
                            BSIM3cgdb = 0.0;

                            BSIM3cdgb = 0.0;
                            BSIM3cddb = 0.0;
                            BSIM3cdsb = 0.0;

                            BSIM3cbgb = -(BSIM3cggb
                                                - Two_Third_CoxWL * dVgs_eff_dVg);
                            T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
                            BSIM3cbsb = -(BSIM3cbgb + T3);
                            BSIM3cbdb = 0.0;
                            BSIM3qinv = -(qgate + qbulk);
                        }
                        else
                        {   /* linear region */
                            Alphaz = Vgst / Vdsat;
                            T1 = 2.0 * Vdsat - Vds;
                            T2 = Vds / (3.0 * T1);
                            T3 = T2 * Vds;
                            T9 = 0.25 * CoxWL;
                            T4 = T9 * Alphaz;
                            T7 = 2.0 * Vds - T1 - 3.0 * T3;
                            T8 = T3 - T1 - 2.0 * Vds;
                            qgate = CoxWL * (Vgs_eff - Vfb
                                             - param.BSIM3phi - 0.5 * (Vds - T3));
                            T10 = T4 * T8;
                            qdrn = T4 * T7;
                            qbulk = -(qgate + qdrn + T10);

                            T5 = T3 / T1;
                            BSIM3cggb = CoxWL * (1.0 - T5 * dVdsat_dVg)
                                              * dVgs_eff_dVg;
                            T11 = -CoxWL * T5 * dVdsat_dVb;
                            BSIM3cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
                            BSIM3cgsb = -(BSIM3cggb + T11
                                                + BSIM3cgdb);
                            T6 = 1.0 / Vdsat;
                            dAlphaz_dVg = T6 * (1.0 - Alphaz * dVdsat_dVg);
                            dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);
                            T7 = T9 * T7;
                            T8 = T9 * T8;
                            T9 = 2.0 * T4 * (1.0 - 3.0 * T5);
                            BSIM3cdgb = (T7 * dAlphaz_dVg - T9
                                               * dVdsat_dVg) * dVgs_eff_dVg;
                            T12 = T7 * dAlphaz_dVb - T9 * dVdsat_dVb;
                            BSIM3cddb = T4 * (3.0 - 6.0 * T2 - 3.0 * T5);
                            BSIM3cdsb = -(BSIM3cdgb + T12
                                                + BSIM3cddb);

                            T9 = 2.0 * T4 * (1.0 + T5);
                            T10 = (T8 * dAlphaz_dVg - T9 * dVdsat_dVg)
                                  * dVgs_eff_dVg;
                            T11 = T8 * dAlphaz_dVb - T9 * dVdsat_dVb;
                            T12 = T4 * (2.0 * T2 + T5 - 1.0);
                            T0 = -(T10 + T11 + T12);

                            BSIM3cbgb = -(BSIM3cggb
                                                + BSIM3cdgb + T10);
                            BSIM3cbdb = -(BSIM3cgdb
                                                + BSIM3cddb + T12);
                            BSIM3cbsb = -(BSIM3cgsb
                                                + BSIM3cdsb + T0);
                            BSIM3qinv = -(qgate + qbulk);
                        }
                    }
                    else if (model.BSIM3xpart < 0.5)
                    {   /* 40/60 Charge partition model */
                        if (Vds >= Vdsat)
                        {   /* saturation region */
                            T1 = Vdsat / 3.0;
                            qgate = CoxWL * (Vgs_eff - Vfb
                                             - param.BSIM3phi - T1);
                            T2 = -Two_Third_CoxWL * Vgst;
                            qbulk = -(qgate + T2);
                            qdrn = 0.4 * T2;

                            BSIM3cggb = One_Third_CoxWL * (3.0
                                                                 - dVdsat_dVg) * dVgs_eff_dVg;
                            T2 = -One_Third_CoxWL * dVdsat_dVb;
                            BSIM3cgsb = -(BSIM3cggb + T2);
                            BSIM3cgdb = 0.0;

                            T3 = 0.4 * Two_Third_CoxWL;
                            BSIM3cdgb = -T3 * dVgs_eff_dVg;
                            BSIM3cddb = 0.0;
                            T4 = T3 * dVth_dVb;
                            BSIM3cdsb = -(T4 + BSIM3cdgb);

                            BSIM3cbgb = -(BSIM3cggb
                                                - Two_Third_CoxWL * dVgs_eff_dVg);
                            T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
                            BSIM3cbsb = -(BSIM3cbgb + T3);
                            BSIM3cbdb = 0.0;
                            BSIM3qinv = -(qgate + qbulk);
                        }
                        else
                        {   /* linear region  */
                            Alphaz = Vgst / Vdsat;
                            T1 = 2.0 * Vdsat - Vds;
                            T2 = Vds / (3.0 * T1);
                            T3 = T2 * Vds;
                            T9 = 0.25 * CoxWL;
                            T4 = T9 * Alphaz;
                            qgate = CoxWL * (Vgs_eff - Vfb - param.BSIM3phi
                                             - 0.5 * (Vds - T3));

                            T5 = T3 / T1;
                            BSIM3cggb = CoxWL * (1.0 - T5 * dVdsat_dVg)
                                              * dVgs_eff_dVg;
                            tmp = -CoxWL * T5 * dVdsat_dVb;
                            BSIM3cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
                            BSIM3cgsb = -(BSIM3cggb
                                                + BSIM3cgdb + tmp);

                            T6 = 1.0 / Vdsat;
                            dAlphaz_dVg = T6 * (1.0 - Alphaz * dVdsat_dVg);
                            dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);

                            T6 = 8.0 * Vdsat * Vdsat - 6.0 * Vdsat * Vds
                                 + 1.2 * Vds * Vds;
                            T8 = T2 / T1;
                            T7 = Vds - T1 - T8 * T6;
                            qdrn = T4 * T7;
                            T7 *= T9;
                            tmp = T8 / T1;
                            tmp1 = T4 * (2.0 - 4.0 * tmp * T6
                                         + T8 * (16.0 * Vdsat - 6.0 * Vds));

                            BSIM3cdgb = (T7 * dAlphaz_dVg - tmp1
                                               * dVdsat_dVg) * dVgs_eff_dVg;
                            T10 = T7 * dAlphaz_dVb - tmp1 * dVdsat_dVb;
                            BSIM3cddb = T4 * (2.0 - (1.0 / (3.0 * T1
                                                                  * T1) + 2.0 * tmp) * T6 + T8
                                                    * (6.0 * Vdsat - 2.4 * Vds));
                            BSIM3cdsb = -(BSIM3cdgb
                                                + T10 + BSIM3cddb);

                            T7 = 2.0 * (T1 + T3);
                            qbulk = -(qgate - T4 * T7);
                            T7 *= T9;
                            T0 = 4.0 * T4 * (1.0 - T5);
                            T12 = (-T7 * dAlphaz_dVg - BSIM3cdgb
                                   - T0 * dVdsat_dVg) * dVgs_eff_dVg;
                            T11 = -T7 * dAlphaz_dVb - T10 - T0 * dVdsat_dVb;
                            T10 = -4.0 * T4 * (T2 - 0.5 + 0.5 * T5)
                                  - BSIM3cddb;
                            tmp = -(T10 + T11 + T12);

                            BSIM3cbgb = -(BSIM3cggb
                                                + BSIM3cdgb + T12);
                            BSIM3cbdb = -(BSIM3cgdb
                                                + BSIM3cddb + T10); /* bug fix */
                            BSIM3cbsb = -(BSIM3cgsb
                                                + BSIM3cdsb + tmp);
                            BSIM3qinv = -(qgate + qbulk);
                        }
                    }
                    else
                    {   /* 50/50 partitioning */
                        if (Vds >= Vdsat)
                        {   /* saturation region */
                            T1 = Vdsat / 3.0;
                            qgate = CoxWL * (Vgs_eff - Vfb
                                             - param.BSIM3phi - T1);
                            T2 = -Two_Third_CoxWL * Vgst;
                            qbulk = -(qgate + T2);
                            qdrn = 0.5 * T2;

                            BSIM3cggb = One_Third_CoxWL * (3.0
                                                                 - dVdsat_dVg) * dVgs_eff_dVg;
                            T2 = -One_Third_CoxWL * dVdsat_dVb;
                            BSIM3cgsb = -(BSIM3cggb + T2);
                            BSIM3cgdb = 0.0;

                            BSIM3cdgb = -One_Third_CoxWL * dVgs_eff_dVg;
                            BSIM3cddb = 0.0;
                            T4 = One_Third_CoxWL * dVth_dVb;
                            BSIM3cdsb = -(T4 + BSIM3cdgb);

                            BSIM3cbgb = -(BSIM3cggb
                                                - Two_Third_CoxWL * dVgs_eff_dVg);
                            T3 = -(T2 + Two_Third_CoxWL * dVth_dVb);
                            BSIM3cbsb = -(BSIM3cbgb + T3);
                            BSIM3cbdb = 0.0;
                            BSIM3qinv = -(qgate + qbulk);
                        }
                        else
                        {   /* linear region */
                            Alphaz = Vgst / Vdsat;
                            T1 = 2.0 * Vdsat - Vds;
                            T2 = Vds / (3.0 * T1);
                            T3 = T2 * Vds;
                            T9 = 0.25 * CoxWL;
                            T4 = T9 * Alphaz;
                            qgate = CoxWL * (Vgs_eff - Vfb - param.BSIM3phi
                                             - 0.5 * (Vds - T3));

                            T5 = T3 / T1;
                            BSIM3cggb = CoxWL * (1.0 - T5 * dVdsat_dVg)
                                              * dVgs_eff_dVg;
                            tmp = -CoxWL * T5 * dVdsat_dVb;
                            BSIM3cgdb = CoxWL * (T2 - 0.5 + 0.5 * T5);
                            BSIM3cgsb = -(BSIM3cggb
                                                + BSIM3cgdb + tmp);

                            T6 = 1.0 / Vdsat;
                            dAlphaz_dVg = T6 * (1.0 - Alphaz * dVdsat_dVg);
                            dAlphaz_dVb = -T6 * (dVth_dVb + Alphaz * dVdsat_dVb);

                            T7 = T1 + T3;
                            qdrn = -T4 * T7;
                            qbulk = -(qgate + qdrn + qdrn);
                            T7 *= T9;
                            T0 = T4 * (2.0 * T5 - 2.0);

                            BSIM3cdgb = (T0 * dVdsat_dVg - T7
                                               * dAlphaz_dVg) * dVgs_eff_dVg;
                            T12 = T0 * dVdsat_dVb - T7 * dAlphaz_dVb;
                            BSIM3cddb = T4 * (1.0 - 2.0 * T2 - T5);
                            BSIM3cdsb = -(BSIM3cdgb + T12
                                                + BSIM3cddb);

                            BSIM3cbgb = -(BSIM3cggb
                                                + 2.0 * BSIM3cdgb);
                            BSIM3cbdb = -(BSIM3cgdb
                                                + 2.0 * BSIM3cddb);
                            BSIM3cbsb = -(BSIM3cgsb
                                                + 2.0 * BSIM3cdsb);
                            BSIM3qinv = -(qgate + qbulk);
                        }
                    }
                }
            }
            else
            {
                if (Vbseff < 0.0)
                {
                    VbseffCV = Vbseff;
                    dVbseffCV_dVb = 1.0;
                }
                else
                {
                    VbseffCV = param.BSIM3phi - Phis;
                    dVbseffCV_dVb = -dPhis_dVb;
                }

                CoxWL = model.BSIM3cox * param.BSIM3weffCV
                        * param.BSIM3leffCV;

                /* Seperate VgsteffCV with noff and voffcv */
                noff = n * param.BSIM3noff;
                dnoff_dVd = param.BSIM3noff * dn_dVd;
                dnoff_dVb = param.BSIM3noff * dn_dVb;
                T0 = Vtm * noff;
                voffcv = param.BSIM3voffcv;
                VgstNVt = (Vgst - voffcv) / T0;

                if (VgstNVt > BSIM3Model.EXP_THRESHOLD)
                {
                    Vgsteff = Vgst - voffcv;
                    dVgsteff_dVg = dVgs_eff_dVg;
                    dVgsteff_dVd = -dVth_dVd;
                    dVgsteff_dVb = -dVth_dVb;
                }
                else if (VgstNVt < -BSIM3Model.EXP_THRESHOLD)
                {
                    Vgsteff = T0 * Math.Log(1.0 + BSIM3Model.MIN_EXP);
                    dVgsteff_dVg = 0.0;
                    dVgsteff_dVd = Vgsteff / noff;
                    dVgsteff_dVb = dVgsteff_dVd * dnoff_dVb;
                    dVgsteff_dVd *= dnoff_dVd;
                }
                else
                {
                    ExpVgst = Math.Exp(VgstNVt);
                    Vgsteff = T0 * Math.Log(1.0 + ExpVgst);
                    dVgsteff_dVg = ExpVgst / (1.0 + ExpVgst);
                    dVgsteff_dVd = -dVgsteff_dVg * (dVth_dVd + (Vgst - voffcv)
                                                    / noff * dnoff_dVd) + Vgsteff / noff * dnoff_dVd;
                    dVgsteff_dVb = -dVgsteff_dVg * (dVth_dVb + (Vgst - voffcv)
                                                    / noff * dnoff_dVb) + Vgsteff / noff * dnoff_dVb;
                    dVgsteff_dVg *= dVgs_eff_dVg;
                } /* End of VgsteffCV */

                if (model.BSIM3capMod == 1)
                {
                    Vfb = param.BSIM3vfbzb;
                    Arg1 = Vgs_eff - VbseffCV - Vfb - Vgsteff;

                    if (Arg1 <= 0.0)
                    {
                        qgate = CoxWL * Arg1;
                        Cgg = CoxWL * (dVgs_eff_dVg - dVgsteff_dVg);
                        Cgd = -CoxWL * dVgsteff_dVd;
                        Cgb = -CoxWL * (dVbseffCV_dVb + dVgsteff_dVb);
                    }
                    else
                    {
                        T0 = 0.5 * param.BSIM3k1ox;
                        T1 = Math.Sqrt(T0 * T0 + Arg1);
                        T2 = CoxWL * T0 / T1;

                        qgate = CoxWL * param.BSIM3k1ox * (T1 - T0);

                        Cgg = T2 * (dVgs_eff_dVg - dVgsteff_dVg);
                        Cgd = -T2 * dVgsteff_dVd;
                        Cgb = -T2 * (dVbseffCV_dVb + dVgsteff_dVb);
                    }
                    qbulk = -qgate;
                    Cbg = -Cgg;
                    Cbd = -Cgd;
                    Cbb = -Cgb;

                    One_Third_CoxWL = CoxWL / 3.0;
                    Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
                    AbulkCV = Abulk0 * param.BSIM3abulkCVfactor;
                    dAbulkCV_dVb = param.BSIM3abulkCVfactor * dAbulk0_dVb;
                    VdsatCV = Vgsteff / AbulkCV;
                    if (VdsatCV < Vds)
                    {
                        dVdsatCV_dVg = 1.0 / AbulkCV;
                        dVdsatCV_dVb = -VdsatCV * dAbulkCV_dVb / AbulkCV;
                        T0 = Vgsteff - VdsatCV / 3.0;
                        dT0_dVg = 1.0 - dVdsatCV_dVg / 3.0;
                        dT0_dVb = -dVdsatCV_dVb / 3.0;
                        qgate += CoxWL * T0;
                        Cgg1 = CoxWL * dT0_dVg;
                        Cgb1 = CoxWL * dT0_dVb + Cgg1 * dVgsteff_dVb;
                        Cgd1 = Cgg1 * dVgsteff_dVd;
                        Cgg1 *= dVgsteff_dVg;
                        Cgg += Cgg1;
                        Cgb += Cgb1;
                        Cgd += Cgd1;

                        T0 = VdsatCV - Vgsteff;
                        dT0_dVg = dVdsatCV_dVg - 1.0;
                        dT0_dVb = dVdsatCV_dVb;
                        qbulk += One_Third_CoxWL * T0;
                        Cbg1 = One_Third_CoxWL * dT0_dVg;
                        Cbb1 = One_Third_CoxWL * dT0_dVb + Cbg1 * dVgsteff_dVb;
                        Cbd1 = Cbg1 * dVgsteff_dVd;
                        Cbg1 *= dVgsteff_dVg;
                        Cbg += Cbg1;
                        Cbb += Cbb1;
                        Cbd += Cbd1;

                        if (model.BSIM3xpart > 0.5)
                            T0 = -Two_Third_CoxWL;
                        else if (model.BSIM3xpart < 0.5)
                            T0 = -0.4 * CoxWL;
                        else
                            T0 = -One_Third_CoxWL;

                        qsrc = T0 * Vgsteff;
                        Csg = T0 * dVgsteff_dVg;
                        Csb = T0 * dVgsteff_dVb;
                        Csd = T0 * dVgsteff_dVd;
                        Cgb *= dVbseff_dVb;
                        Cbb *= dVbseff_dVb;
                        Csb *= dVbseff_dVb;
                    }
                    else
                    {
                        T0 = AbulkCV * Vds;
                        T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1.0e-20);
                        T2 = Vds / T1;
                        T3 = T0 * T2;
                        dT3_dVg = -12.0 * T2 * T2 * AbulkCV;
                        dT3_dVd = 6.0 * T0 * (4.0 * Vgsteff - T0) / T1 / T1 - 0.5;
                        dT3_dVb = 12.0 * T2 * T2 * dAbulkCV_dVb * Vgsteff;

                        qgate += CoxWL * (Vgsteff - 0.5 * Vds + T3);
                        Cgg1 = CoxWL * (1.0 + dT3_dVg);
                        Cgb1 = CoxWL * dT3_dVb + Cgg1 * dVgsteff_dVb;
                        Cgd1 = CoxWL * dT3_dVd + Cgg1 * dVgsteff_dVd;
                        Cgg1 *= dVgsteff_dVg;
                        Cgg += Cgg1;
                        Cgb += Cgb1;
                        Cgd += Cgd1;

                        qbulk += CoxWL * (1.0 - AbulkCV) * (0.5 * Vds - T3);
                        Cbg1 = -CoxWL * ((1.0 - AbulkCV) * dT3_dVg);
                        Cbb1 = -CoxWL * ((1.0 - AbulkCV) * dT3_dVb
                                         + (0.5 * Vds - T3) * dAbulkCV_dVb)
                               + Cbg1 * dVgsteff_dVb;
                        Cbd1 = -CoxWL * (1.0 - AbulkCV) * dT3_dVd
                               + Cbg1 * dVgsteff_dVd;
                        Cbg1 *= dVgsteff_dVg;
                        Cbg += Cbg1;
                        Cbb += Cbb1;
                        Cbd += Cbd1;

                        if (model.BSIM3xpart > 0.5)
                        {   /* 0/100 Charge petition model */
                            T1 = T1 + T1;
                            qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0
                                             - T0 * T0 / T1);
                            Csg = -CoxWL * (0.5 + 24.0 * T0 * Vds / T1 / T1
                                            * AbulkCV);
                            Csb = -CoxWL * (0.25 * Vds * dAbulkCV_dVb
                                            - 12.0 * T0 * Vds / T1 / T1 * (4.0 * Vgsteff - T0)
                                            * dAbulkCV_dVb) + Csg * dVgsteff_dVb;
                            Csd = -CoxWL * (0.25 * AbulkCV - 12.0 * AbulkCV * T0
                                            / T1 / T1 * (4.0 * Vgsteff - T0))
                                  + Csg * dVgsteff_dVd;
                            Csg *= dVgsteff_dVg;
                        }
                        else if (model.BSIM3xpart < 0.5)
                        {   /* 40/60 Charge petition model */
                            T1 = T1 / 12.0;
                            T2 = 0.5 * CoxWL / (T1 * T1);
                            T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff
                                            * (Vgsteff - 4.0 * T0 / 3.0))
                                 - 2.0 * T0 * T0 * T0 / 15.0;
                            qsrc = -T2 * T3;
                            T4 = 4.0 / 3.0 * Vgsteff * (Vgsteff - T0)
                                 + 0.4 * T0 * T0;
                            Csg = -2.0 * qsrc / T1 - T2 * (Vgsteff * (3.0
                                                           * Vgsteff - 8.0 * T0 / 3.0)
                                                           + 2.0 * T0 * T0 / 3.0);
                            Csb = (qsrc / T1 * Vds + T2 * T4 * Vds) * dAbulkCV_dVb
                                  + Csg * dVgsteff_dVb;
                            Csd = (qsrc / T1 + T2 * T4) * AbulkCV
                                  + Csg * dVgsteff_dVd;
                            Csg *= dVgsteff_dVg;
                        }
                        else
                        {   /* 50/50 Charge petition model */
                            qsrc = -0.5 * (qgate + qbulk);
                            Csg = -0.5 * (Cgg1 + Cbg1);
                            Csb = -0.5 * (Cgb1 + Cbb1);
                            Csd = -0.5 * (Cgd1 + Cbd1);
                        }
                        Cgb *= dVbseff_dVb;
                        Cbb *= dVbseff_dVb;
                        Csb *= dVbseff_dVb;
                    }
                    qdrn = -(qgate + qbulk + qsrc);
                    BSIM3cggb = Cgg;
                    BSIM3cgsb = -(Cgg + Cgd + Cgb);
                    BSIM3cgdb = Cgd;
                    BSIM3cdgb = -(Cgg + Cbg + Csg);
                    BSIM3cdsb = (Cgg + Cgd + Cgb + Cbg + Cbd + Cbb
                                       + Csg + Csd + Csb);
                    BSIM3cddb = -(Cgd + Cbd + Csd);
                    BSIM3cbgb = Cbg;
                    BSIM3cbsb = -(Cbg + Cbd + Cbb);
                    BSIM3cbdb = Cbd;
                    BSIM3qinv = -(qgate + qbulk);
                }

                else if (model.BSIM3capMod == 2)
                {
                    Vfb = param.BSIM3vfbzb;
                    V3 = Vfb - Vgs_eff + VbseffCV - DELTA_3;
                    if (Vfb <= 0.0)
                    {
                        T0 = Math.Sqrt(V3 * V3 - 4.0 * DELTA_3 * Vfb);
                        T2 = -DELTA_3 / T0;
                    }
                    else
                    {
                        T0 = Math.Sqrt(V3 * V3 + 4.0 * DELTA_3 * Vfb);
                        T2 = DELTA_3 / T0;
                    }

                    T1 = 0.5 * (1.0 + V3 / T0);
                    Vfbeff = Vfb - 0.5 * (V3 + T0);
                    dVfbeff_dVg = T1 * dVgs_eff_dVg;
                    dVfbeff_dVb = -T1 * dVbseffCV_dVb;
                    Qac0 = CoxWL * (Vfbeff - Vfb);
                    dQac0_dVg = CoxWL * dVfbeff_dVg;
                    dQac0_dVb = CoxWL * dVfbeff_dVb;

                    T0 = 0.5 * param.BSIM3k1ox;
                    T3 = Vgs_eff - Vfbeff - VbseffCV - Vgsteff;
                    if (param.BSIM3k1ox == 0.0)
                    {
                        T1 = 0.0;
                        T2 = 0.0;
                    }
                    else if (T3 < 0.0)
                    {
                        T1 = T0 + T3 / param.BSIM3k1ox;
                        T2 = CoxWL;
                    }
                    else
                    {
                        T1 = Math.Sqrt(T0 * T0 + T3);
                        T2 = CoxWL * T0 / T1;
                    }

                    Qsub0 = CoxWL * param.BSIM3k1ox * (T1 - T0);

                    dQsub0_dVg = T2 * (dVgs_eff_dVg - dVfbeff_dVg - dVgsteff_dVg);
                    dQsub0_dVd = -T2 * dVgsteff_dVd;
                    dQsub0_dVb = -T2 * (dVfbeff_dVb + dVbseffCV_dVb
                                        + dVgsteff_dVb);

                    AbulkCV = Abulk0 * param.BSIM3abulkCVfactor;
                    dAbulkCV_dVb = param.BSIM3abulkCVfactor * dAbulk0_dVb;
                    VdsatCV = Vgsteff / AbulkCV;

                    V4 = VdsatCV - Vds - DELTA_4;
                    T0 = Math.Sqrt(V4 * V4 + 4.0 * DELTA_4 * VdsatCV);
                    VdseffCV = VdsatCV - 0.5 * (V4 + T0);
                    T1 = 0.5 * (1.0 + V4 / T0);
                    T2 = DELTA_4 / T0;
                    T3 = (1.0 - T1 - T2) / AbulkCV;
                    dVdseffCV_dVg = T3;
                    dVdseffCV_dVd = T1;
                    dVdseffCV_dVb = -T3 * VdsatCV * dAbulkCV_dVb;
                    /* Added to eliminate non-zero VdseffCV at Vds=0.0 */
                    if (Vds == 0.0)
                    {
                        VdseffCV = 0.0;
                        dVdseffCV_dVg = 0.0;
                        dVdseffCV_dVb = 0.0;
                    }

                    T0 = AbulkCV * VdseffCV;
                    T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1e-20);
                    T2 = VdseffCV / T1;
                    T3 = T0 * T2;

                    T4 = (1.0 - 12.0 * T2 * T2 * AbulkCV);
                    T5 = (6.0 * T0 * (4.0 * Vgsteff - T0) / (T1 * T1) - 0.5);
                    T6 = 12.0 * T2 * T2 * Vgsteff;

                    qinoi = -CoxWL * (Vgsteff - 0.5 * T0 + AbulkCV * T3);
                    qgate = CoxWL * (Vgsteff - 0.5 * VdseffCV + T3);
                    Cgg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
                    Cgd1 = CoxWL * T5 * dVdseffCV_dVd + Cgg1 * dVgsteff_dVd;
                    Cgb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                           + Cgg1 * dVgsteff_dVb;
                    Cgg1 *= dVgsteff_dVg;

                    T7 = 1.0 - AbulkCV;
                    qbulk = CoxWL * T7 * (0.5 * VdseffCV - T3);
                    T4 = -T7 * (T4 - 1.0);
                    T5 = -T7 * T5;
                    T6 = -(T7 * T6 + (0.5 * VdseffCV - T3));
                    Cbg1 = CoxWL * (T4 + T5 * dVdseffCV_dVg);
                    Cbd1 = CoxWL * T5 * dVdseffCV_dVd + Cbg1 * dVgsteff_dVd;
                    Cbb1 = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                           + Cbg1 * dVgsteff_dVb;
                    Cbg1 *= dVgsteff_dVg;

                    if (model.BSIM3xpart > 0.5)
                    {   /* 0/100 Charge petition model */
                        T1 = T1 + T1;
                        qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0
                                         - T0 * T0 / T1);
                        T7 = (4.0 * Vgsteff - T0) / (T1 * T1);
                        T4 = -(0.5 + 24.0 * T0 * T0 / (T1 * T1));
                        T5 = -(0.25 * AbulkCV - 12.0 * AbulkCV * T0 * T7);
                        T6 = -(0.25 * VdseffCV - 12.0 * T0 * VdseffCV * T7);
                        Csg = CoxWL * (T4 + T5 * dVdseffCV_dVg);
                        Csd = CoxWL * T5 * dVdseffCV_dVd + Csg * dVgsteff_dVd;
                        Csb = CoxWL * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                              + Csg * dVgsteff_dVb;
                        Csg *= dVgsteff_dVg;
                    }
                    else if (model.BSIM3xpart < 0.5)
                    {   /* 40/60 Charge petition model */
                        T1 = T1 / 12.0;
                        T2 = 0.5 * CoxWL / (T1 * T1);
                        T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff
                                        * (Vgsteff - 4.0 * T0 / 3.0))
                             - 2.0 * T0 * T0 * T0 / 15.0;
                        qsrc = -T2 * T3;
                        T7 = 4.0 / 3.0 * Vgsteff * (Vgsteff - T0)
                             + 0.4 * T0 * T0;
                        T4 = -2.0 * qsrc / T1 - T2 * (Vgsteff * (3.0
                                                      * Vgsteff - 8.0 * T0 / 3.0)
                                                      + 2.0 * T0 * T0 / 3.0);
                        T5 = (qsrc / T1 + T2 * T7) * AbulkCV;
                        T6 = (qsrc / T1 * VdseffCV + T2 * T7 * VdseffCV);
                        Csg = (T4 + T5 * dVdseffCV_dVg);
                        Csd = T5 * dVdseffCV_dVd + Csg * dVgsteff_dVd;
                        Csb = (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                              + Csg * dVgsteff_dVb;
                        Csg *= dVgsteff_dVg;
                    }
                    else
                    {   /* 50/50 Charge petition model */
                        qsrc = -0.5 * (qgate + qbulk);
                        Csg = -0.5 * (Cgg1 + Cbg1);
                        Csb = -0.5 * (Cgb1 + Cbb1);
                        Csd = -0.5 * (Cgd1 + Cbd1);
                    }

                    qgate += Qac0 + Qsub0;
                    qbulk -= (Qac0 + Qsub0);
                    qdrn = -(qgate + qbulk + qsrc);

                    Cgg = dQac0_dVg + dQsub0_dVg + Cgg1;
                    Cgd = dQsub0_dVd + Cgd1;
                    Cgb = dQac0_dVb + dQsub0_dVb + Cgb1;

                    Cbg = Cbg1 - dQac0_dVg - dQsub0_dVg;
                    Cbd = Cbd1 - dQsub0_dVd;
                    Cbb = Cbb1 - dQac0_dVb - dQsub0_dVb;

                    Cgb *= dVbseff_dVb;
                    Cbb *= dVbseff_dVb;
                    Csb *= dVbseff_dVb;

                    BSIM3cggb = Cgg;
                    BSIM3cgsb = -(Cgg + Cgd + Cgb);
                    BSIM3cgdb = Cgd;
                    BSIM3cdgb = -(Cgg + Cbg + Csg);
                    BSIM3cdsb = (Cgg + Cgd + Cgb + Cbg + Cbd + Cbb
                                       + Csg + Csd + Csb);
                    BSIM3cddb = -(Cgd + Cbd + Csd);
                    BSIM3cbgb = Cbg;
                    BSIM3cbsb = -(Cbg + Cbd + Cbb);
                    BSIM3cbdb = Cbd;
                    BSIM3qinv = qinoi;
                }

                /* New Charge-Thickness capMod (CTM) begins */
                else if (model.BSIM3capMod == 3)
                {
                    V3 = param.BSIM3vfbzb - Vgs_eff + VbseffCV - DELTA_3;
                    if (param.BSIM3vfbzb <= 0.0)
                    {
                        T0 = Math.Sqrt(V3 * V3 - 4.0 * DELTA_3 * param.BSIM3vfbzb);
                        T2 = -DELTA_3 / T0;
                    }
                    else
                    {
                        T0 = Math.Sqrt(V3 * V3 + 4.0 * DELTA_3 * param.BSIM3vfbzb);
                        T2 = DELTA_3 / T0;
                    }

                    T1 = 0.5 * (1.0 + V3 / T0);
                    Vfbeff = param.BSIM3vfbzb - 0.5 * (V3 + T0);
                    dVfbeff_dVg = T1 * dVgs_eff_dVg;
                    dVfbeff_dVb = -T1 * dVbseffCV_dVb;

                    Cox = model.BSIM3cox;
                    Tox = 1.0e8 * model.BSIM3tox;
                    T0 = (Vgs_eff - VbseffCV - param.BSIM3vfbzb) / Tox;
                    dT0_dVg = dVgs_eff_dVg / Tox;
                    dT0_dVb = -dVbseffCV_dVb / Tox;

                    tmp = T0 * param.BSIM3acde;
                    if ((-BSIM3Model.EXP_THRESHOLD < tmp) && (tmp < BSIM3Model.EXP_THRESHOLD))
                    {
                        Tcen = param.BSIM3ldeb * Math.Exp(tmp);
                        dTcen_dVg = param.BSIM3acde * Tcen;
                        dTcen_dVb = dTcen_dVg * dT0_dVb;
                        dTcen_dVg *= dT0_dVg;
                    }
                    else if (tmp <= -BSIM3Model.EXP_THRESHOLD)
                    {
                        Tcen = param.BSIM3ldeb * BSIM3Model.MIN_EXP;
                        dTcen_dVg = dTcen_dVb = 0.0;
                    }
                    else
                    {
                        Tcen = param.BSIM3ldeb * BSIM3Model.MAX_EXP;
                        dTcen_dVg = dTcen_dVb = 0.0;
                    }

                    LINK = 1.0e-3 * model.BSIM3tox;
                    V3 = param.BSIM3ldeb - Tcen - LINK;
                    V4 = Math.Sqrt(V3 * V3 + 4.0 * LINK * param.BSIM3ldeb);
                    Tcen = param.BSIM3ldeb - 0.5 * (V3 + V4);
                    T1 = 0.5 * (1.0 + V3 / V4);
                    dTcen_dVg *= T1;
                    dTcen_dVb *= T1;

                    Ccen = BSIM3Model.EPSSI / Tcen;
                    T2 = Cox / (Cox + Ccen);
                    Coxeff = T2 * Ccen;
                    T3 = -Ccen / Tcen;
                    dCoxeff_dVg = T2 * T2 * T3;
                    dCoxeff_dVb = dCoxeff_dVg * dTcen_dVb;
                    dCoxeff_dVg *= dTcen_dVg;
                    CoxWLcen = CoxWL * Coxeff / Cox;

                    Qac0 = CoxWLcen * (Vfbeff - param.BSIM3vfbzb);
                    QovCox = Qac0 / Coxeff;
                    dQac0_dVg = CoxWLcen * dVfbeff_dVg
                                + QovCox * dCoxeff_dVg;
                    dQac0_dVb = CoxWLcen * dVfbeff_dVb
                                + QovCox * dCoxeff_dVb;

                    T0 = 0.5 * param.BSIM3k1ox;
                    T3 = Vgs_eff - Vfbeff - VbseffCV - Vgsteff;
                    if (param.BSIM3k1ox == 0.0)
                    {
                        T1 = 0.0;
                        T2 = 0.0;
                    }
                    else if (T3 < 0.0)
                    {
                        T1 = T0 + T3 / param.BSIM3k1ox;
                        T2 = CoxWLcen;
                    }
                    else
                    {
                        T1 = Math.Sqrt(T0 * T0 + T3);
                        T2 = CoxWLcen * T0 / T1;
                    }

                    Qsub0 = CoxWLcen * param.BSIM3k1ox * (T1 - T0);
                    QovCox = Qsub0 / Coxeff;
                    dQsub0_dVg = T2 * (dVgs_eff_dVg - dVfbeff_dVg - dVgsteff_dVg)
                                 + QovCox * dCoxeff_dVg;
                    dQsub0_dVd = -T2 * dVgsteff_dVd;
                    dQsub0_dVb = -T2 * (dVfbeff_dVb + dVbseffCV_dVb + dVgsteff_dVb)
                                 + QovCox * dCoxeff_dVb;

                    /* Gate-bias dependent delta Phis begins */
                    if (param.BSIM3k1ox <= 0.0)
                    {
                        Denomi = 0.25 * param.BSIM3moin * Vtm;
                        T0 = 0.5 * param.BSIM3sqrtPhi;
                    }
                    else
                    {
                        Denomi = param.BSIM3moin * Vtm
                                 * param.BSIM3k1ox * param.BSIM3k1ox;
                        T0 = param.BSIM3k1ox * param.BSIM3sqrtPhi;
                    }
                    T1 = 2.0 * T0 + Vgsteff;

                    DeltaPhi = Vtm * Math.Log(1.0 + T1 * Vgsteff / Denomi);
                    dDeltaPhi_dVg = 2.0 * Vtm * (T1 - T0) / (Denomi + T1 * Vgsteff);
                    /* End of delta Phis */

                    /* VgDP = Vgsteff - DeltaPhi */
                    T0 = Vgsteff - DeltaPhi - 0.001;
                    dT0_dVg = 1.0 - dDeltaPhi_dVg;
                    T1 = Math.Sqrt(T0 * T0 + Vgsteff * 0.004);
                    VgDP = 0.5 * (T0 + T1);
                    dVgDP_dVg = 0.5 * (dT0_dVg + (T0 * dT0_dVg + 0.002) / T1);

                    T3 = 4.0 * (Vth - param.BSIM3vfbzb - param.BSIM3phi);
                    Tox += Tox;
                    if (T3 >= 0.0)
                    {
                        T0 = (Vgsteff + T3) / Tox;
                        dT0_dVd = (dVgsteff_dVd + 4.0 * dVth_dVd) / Tox;
                        dT0_dVb = (dVgsteff_dVb + 4.0 * dVth_dVb) / Tox;
                    }
                    else
                    {
                        T0 = (Vgsteff + 1.0e-20) / Tox;
                        dT0_dVd = dVgsteff_dVd / Tox;
                        dT0_dVb = dVgsteff_dVb / Tox;
                    }
                    tmp = Math.Exp(0.7 * Math.Log(T0));
                    T1 = 1.0 + tmp;
                    T2 = 0.7 * tmp / (T0 * Tox);
                    Tcen = 1.9e-9 / T1;
                    dTcen_dVg = -1.9e-9 * T2 / T1 / T1;
                    dTcen_dVd = Tox * dTcen_dVg;
                    dTcen_dVb = dTcen_dVd * dT0_dVb;
                    dTcen_dVd *= dT0_dVd;
                    dTcen_dVg *= dVgsteff_dVg;

                    Ccen = BSIM3Model.EPSSI / Tcen;
                    T0 = Cox / (Cox + Ccen);
                    Coxeff = T0 * Ccen;
                    T1 = -Ccen / Tcen;
                    dCoxeff_dVg = T0 * T0 * T1;
                    dCoxeff_dVd = dCoxeff_dVg * dTcen_dVd;
                    dCoxeff_dVb = dCoxeff_dVg * dTcen_dVb;
                    dCoxeff_dVg *= dTcen_dVg;
                    CoxWLcen = CoxWL * Coxeff / Cox;

                    AbulkCV = Abulk0 * param.BSIM3abulkCVfactor;
                    dAbulkCV_dVb = param.BSIM3abulkCVfactor * dAbulk0_dVb;
                    VdsatCV = VgDP / AbulkCV;
                    T0 = VdsatCV - Vds - DELTA_4;
                    dT0_dVg = dVgDP_dVg / AbulkCV;
                    dT0_dVb = -VdsatCV * dAbulkCV_dVb / AbulkCV;
                    T1 = Math.Sqrt(T0 * T0 + 4.0 * DELTA_4 * VdsatCV);
                    dT1_dVg = (T0 + DELTA_4 + DELTA_4) / T1;
                    dT1_dVd = -T0 / T1;
                    dT1_dVb = dT1_dVg * dT0_dVb;
                    dT1_dVg *= dT0_dVg;
                    if (T0 >= 0.0)
                    {
                        VdseffCV = VdsatCV - 0.5 * (T0 + T1);
                        dVdseffCV_dVg = 0.5 * (dT0_dVg - dT1_dVg);
                        dVdseffCV_dVd = 0.5 * (1.0 - dT1_dVd);
                        dVdseffCV_dVb = 0.5 * (dT0_dVb - dT1_dVb);
                    }
                    else
                    {
                        T3 = (DELTA_4 + DELTA_4) / (T1 - T0);
                        T4 = 1.0 - T3;
                        T5 = VdsatCV * T3 / (T1 - T0);
                        VdseffCV = VdsatCV * T4;
                        dVdseffCV_dVg = dT0_dVg * T4 + T5 * (dT1_dVg - dT0_dVg);
                        dVdseffCV_dVd = T5 * (dT1_dVd + 1.0);
                        dVdseffCV_dVb = dT0_dVb * (1.0 - T5) + T5 * dT1_dVb;
                    }

                    /* Added to eliminate non-zero VdseffCV at Vds=0.0 */
                    if (Vds == 0.0)
                    {
                        VdseffCV = 0.0;
                        dVdseffCV_dVg = 0.0;
                        dVdseffCV_dVb = 0.0;
                    }

                    T0 = AbulkCV * VdseffCV;
                    T1 = VgDP;
                    T2 = 12.0 * (T1 - 0.5 * T0 + 1.0e-20);
                    T3 = T0 / T2;
                    T4 = 1.0 - 12.0 * T3 * T3;
                    T5 = AbulkCV * (6.0 * T0 * (4.0 * T1 - T0) / (T2 * T2) - 0.5);
                    T6 = T5 * VdseffCV / AbulkCV;

                    qgate = qinoi = CoxWLcen * (T1 - T0 * (0.5 - T3));
                    QovCox = qgate / Coxeff;
                    Cgg1 = CoxWLcen * (T4 * dVgDP_dVg
                                       + T5 * dVdseffCV_dVg);
                    Cgd1 = CoxWLcen * T5 * dVdseffCV_dVd + Cgg1
                           * dVgsteff_dVd + QovCox * dCoxeff_dVd;
                    Cgb1 = CoxWLcen * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                           + Cgg1 * dVgsteff_dVb + QovCox * dCoxeff_dVb;
                    Cgg1 = Cgg1 * dVgsteff_dVg + QovCox * dCoxeff_dVg;


                    T7 = 1.0 - AbulkCV;
                    T8 = T2 * T2;
                    T9 = 12.0 * T7 * T0 * T0 / (T8 * AbulkCV);
                    T10 = T9 * dVgDP_dVg;
                    T11 = -T7 * T5 / AbulkCV;
                    T12 = -(T9 * T1 / AbulkCV + VdseffCV * (0.5 - T0 / T2));

                    qbulk = CoxWLcen * T7 * (0.5 * VdseffCV - T0 * VdseffCV / T2);
                    QovCox = qbulk / Coxeff;
                    Cbg1 = CoxWLcen * (T10 + T11 * dVdseffCV_dVg);
                    Cbd1 = CoxWLcen * T11 * dVdseffCV_dVd + Cbg1
                           * dVgsteff_dVd + QovCox * dCoxeff_dVd;
                    Cbb1 = CoxWLcen * (T11 * dVdseffCV_dVb + T12 * dAbulkCV_dVb)
                           + Cbg1 * dVgsteff_dVb + QovCox * dCoxeff_dVb;
                    Cbg1 = Cbg1 * dVgsteff_dVg + QovCox * dCoxeff_dVg;

                    if (model.BSIM3xpart > 0.5)
                    {   /* 0/100 partition */
                        qsrc = -CoxWLcen * (T1 / 2.0 + T0 / 4.0
                                            - 0.5 * T0 * T0 / T2);
                        QovCox = qsrc / Coxeff;
                        T2 += T2;
                        T3 = T2 * T2;
                        T7 = -(0.25 - 12.0 * T0 * (4.0 * T1 - T0) / T3);
                        T4 = -(0.5 + 24.0 * T0 * T0 / T3) * dVgDP_dVg;
                        T5 = T7 * AbulkCV;
                        T6 = T7 * VdseffCV;

                        Csg = CoxWLcen * (T4 + T5 * dVdseffCV_dVg);
                        Csd = CoxWLcen * T5 * dVdseffCV_dVd + Csg * dVgsteff_dVd
                              + QovCox * dCoxeff_dVd;
                        Csb = CoxWLcen * (T5 * dVdseffCV_dVb + T6 * dAbulkCV_dVb)
                              + Csg * dVgsteff_dVb + QovCox * dCoxeff_dVb;
                        Csg = Csg * dVgsteff_dVg + QovCox * dCoxeff_dVg;
                    }
                    else if (model.BSIM3xpart < 0.5)
                    {   /* 40/60 partition */
                        T2 = T2 / 12.0;
                        T3 = 0.5 * CoxWLcen / (T2 * T2);
                        T4 = T1 * (2.0 * T0 * T0 / 3.0 + T1 * (T1 - 4.0
                                                               * T0 / 3.0)) - 2.0 * T0 * T0 * T0 / 15.0;
                        qsrc = -T3 * T4;
                        QovCox = qsrc / Coxeff;
                        T8 = 4.0 / 3.0 * T1 * (T1 - T0) + 0.4 * T0 * T0;
                        T5 = -2.0 * qsrc / T2 - T3 * (T1 * (3.0 * T1 - 8.0
                                                            * T0 / 3.0) + 2.0 * T0 * T0 / 3.0);
                        T6 = AbulkCV * (qsrc / T2 + T3 * T8);
                        T7 = T6 * VdseffCV / AbulkCV;

                        Csg = T5 * dVgDP_dVg + T6 * dVdseffCV_dVg;
                        Csd = Csg * dVgsteff_dVd + T6 * dVdseffCV_dVd
                              + QovCox * dCoxeff_dVd;
                        Csb = Csg * dVgsteff_dVb + T6 * dVdseffCV_dVb
                              + T7 * dAbulkCV_dVb + QovCox * dCoxeff_dVb;
                        Csg = Csg * dVgsteff_dVg + QovCox * dCoxeff_dVg;
                    }
                    else
                    {   /* 50/50 partition */
                        qsrc = -0.5 * qgate;
                        Csg = -0.5 * Cgg1;
                        Csd = -0.5 * Cgd1;
                        Csb = -0.5 * Cgb1;
                    }

                    qgate += Qac0 + Qsub0 - qbulk;
                    qbulk -= (Qac0 + Qsub0);
                    qdrn = -(qgate + qbulk + qsrc);

                    Cbg = Cbg1 - dQac0_dVg - dQsub0_dVg;
                    Cbd = Cbd1 - dQsub0_dVd;
                    Cbb = Cbb1 - dQac0_dVb - dQsub0_dVb;

                    Cgg = Cgg1 - Cbg;
                    Cgd = Cgd1 - Cbd;
                    Cgb = Cgb1 - Cbb;

                    Cgb *= dVbseff_dVb;
                    Cbb *= dVbseff_dVb;
                    Csb *= dVbseff_dVb;

                    BSIM3cggb = Cgg;
                    BSIM3cgsb = -(Cgg + Cgd + Cgb);
                    BSIM3cgdb = Cgd;
                    BSIM3cdgb = -(Cgg + Cbg + Csg);
                    BSIM3cdsb = (Cgg + Cgd + Cgb + Cbg + Cbd + Cbb
                                       + Csg + Csd + Csb);
                    BSIM3cddb = -(Cgd + Cbd + Csd);
                    BSIM3cbgb = Cbg;
                    BSIM3cbsb = -(Cbg + Cbd + Cbb);
                    BSIM3cbdb = Cbd;
                    BSIM3qinv = -qinoi;
                }  /* End of CTM */
            }

        finished:
            /* Returning Values to Calling Routine */
            /*
             *  COMPUTE EQUIVALENT DRAIN CURRENT SOURCE
             */

            BSIM3qgate = qgate;
            BSIM3qbulk = qbulk;
            BSIM3qdrn = qdrn;
            BSIM3cd = cdrain;

            if (lp.ChargeComputationNeeded)
            {   /*  charge storage elements
               *  bulk-drain and bulk-source depletion capacitances
               *  czbd : zero bias drain junction capacitance
               *  czbs : zero bias source junction capacitance
               *  czbdsw: zero bias drain junction sidewall capacitance
			  along field oxide
               *  czbssw: zero bias source junction sidewall capacitance
			  along field oxide
	       *  czbdswg: zero bias drain junction sidewall capacitance
			   along gate side
	       *  czbsswg: zero bias source junction sidewall capacitance
			   along gate side
               */

                czbd = model.BSIM3unitAreaTempJctCap * BSIM3drainArea; /*bug fix */
                czbs = model.BSIM3unitAreaTempJctCap * BSIM3sourceArea;
                if (BSIM3drainPerimeter < param.BSIM3weff)
                {
                    czbdswg = model.BSIM3unitLengthGateSidewallTempJctCap
                              * BSIM3drainPerimeter;
                    czbdsw = 0.0;
                }
                else
                {
                    czbdsw = model.BSIM3unitLengthSidewallTempJctCap
                             * (BSIM3drainPerimeter - param.BSIM3weff);
                    czbdswg = model.BSIM3unitLengthGateSidewallTempJctCap
                              * param.BSIM3weff;
                }
                if (BSIM3sourcePerimeter < param.BSIM3weff)
                {
                    czbssw = 0.0;
                    czbsswg = model.BSIM3unitLengthGateSidewallTempJctCap
                              * BSIM3sourcePerimeter;
                }
                else
                {
                    czbssw = model.BSIM3unitLengthSidewallTempJctCap
                             * (BSIM3sourcePerimeter - param.BSIM3weff);
                    czbsswg = model.BSIM3unitLengthGateSidewallTempJctCap
                              * param.BSIM3weff;
                }

                MJ = model.BSIM3bulkJctBotGradingCoeff;
                MJSW = model.BSIM3bulkJctSideGradingCoeff;
                MJSWG = model.BSIM3bulkJctGateSideGradingCoeff;

                /* Source Bulk Junction */
                if (vbs == 0.0)
                {
                    state.States[0][BSIM3state + BSIM3qbs] = 0.0;
                    BSIM3capbs = czbs + czbssw + czbsswg;
                }
                else if (vbs < 0.0)
                {
                    if (czbs > 0.0)
                    {
                        arg = 1.0 - vbs / model.BSIM3PhiB;
                        if (MJ == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJ * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbs] = model.BSIM3PhiB * czbs
                                                             * (1.0 - arg * sarg) / (1.0 - MJ);
                        BSIM3capbs = czbs * sarg;
                    }
                    else
                    {
                        state.States[0][BSIM3state + BSIM3qbs] = 0.0;
                        BSIM3capbs = 0.0;
                    }
                    if (czbssw > 0.0)
                    {
                        arg = 1.0 - vbs / model.BSIM3PhiBSW;
                        if (MJSW == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJSW * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbs] += model.BSIM3PhiBSW * czbssw
                                                              * (1.0 - arg * sarg) / (1.0 - MJSW);
                        BSIM3capbs += czbssw * sarg;
                    }
                    if (czbsswg > 0.0)
                    {
                        arg = 1.0 - vbs / model.BSIM3PhiBSWG;
                        if (MJSWG == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJSWG * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbs] += model.BSIM3PhiBSWG * czbsswg
                                                              * (1.0 - arg * sarg) / (1.0 - MJSWG);
                        BSIM3capbs += czbsswg * sarg;
                    }

                }
                else
                {
                    T0 = czbs + czbssw + czbsswg;
                    T1 = vbs * (czbs * MJ / model.BSIM3PhiB + czbssw * MJSW
                                / model.BSIM3PhiBSW + czbsswg * MJSWG / model.BSIM3PhiBSWG);
                    state.States[0][BSIM3state + BSIM3qbs] = vbs * (T0 + 0.5 * T1);
                    BSIM3capbs = T0 + T1;
                }

                /* Drain Bulk Junction */
                if (vbd == 0.0)
                {
                    state.States[0][BSIM3state + BSIM3qbd] = 0.0;
                    BSIM3capbd = czbd + czbdsw + czbdswg;
                }
                else if (vbd < 0.0)
                {
                    if (czbd > 0.0)
                    {
                        arg = 1.0 - vbd / model.BSIM3PhiB;
                        if (MJ == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJ * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbd] = model.BSIM3PhiB * czbd
                                                             * (1.0 - arg * sarg) / (1.0 - MJ);
                        BSIM3capbd = czbd * sarg;
                    }
                    else
                    {
                        state.States[0][BSIM3state + BSIM3qbd] = 0.0;
                        BSIM3capbd = 0.0;
                    }
                    if (czbdsw > 0.0)
                    {
                        arg = 1.0 - vbd / model.BSIM3PhiBSW;
                        if (MJSW == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJSW * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbd] += model.BSIM3PhiBSW * czbdsw
                                                              * (1.0 - arg * sarg) / (1.0 - MJSW);
                        BSIM3capbd += czbdsw * sarg;
                    }
                    if (czbdswg > 0.0)
                    {
                        arg = 1.0 - vbd / model.BSIM3PhiBSWG;
                        if (MJSWG == 0.5)
                            sarg = 1.0 / Math.Sqrt(arg);
                        else
                            sarg = Math.Exp(-MJSWG * Math.Log(arg));
                        state.States[0][BSIM3state + BSIM3qbd] += model.BSIM3PhiBSWG * czbdswg
                                                              * (1.0 - arg * sarg) / (1.0 - MJSWG);
                        BSIM3capbd += czbdswg * sarg;
                    }
                }
                else
                {
                    T0 = czbd + czbdsw + czbdswg;
                    T1 = vbd * (czbd * MJ / model.BSIM3PhiB + czbdsw * MJSW
                                / model.BSIM3PhiBSW + czbdswg * MJSWG / model.BSIM3PhiBSWG);
                    state.States[0][BSIM3state + BSIM3qbd] = vbd * (T0 + 0.5 * T1);
                    BSIM3capbd = T0 + T1;
                }
            }

            /*
             *  check convergence
             */
            if ((!BSIM3off) || (!(ckt.Mode.HasFlag(Circuit.Modes.InitFix))))
            {
                if (Check)
                {
                    ckt.IsCon = false;
                }
            }
            state.States[0][BSIM3state + BSIM3vbs] = vbs;
            state.States[0][BSIM3state + BSIM3vbd] = vbd;
            state.States[0][BSIM3state + BSIM3vgs] = vgs;
            state.States[0][BSIM3state + BSIM3vds] = vds;
            state.States[0][BSIM3state + BSIM3qdef] = qdef;

            /* bulk and channel charge plus overlaps */

            if (!lp.ChargeComputationNeeded)
                goto line850;

            // line755:
            /* NQS begins */
            if ((BSIM3nqsMod != 0) || (BSIM3acnqsMod != 0))
            {
                qcheq = -(qbulk + qgate);

                BSIM3cqgb = -(BSIM3cggb + BSIM3cbgb);
                BSIM3cqdb = -(BSIM3cgdb + BSIM3cbdb);
                BSIM3cqsb = -(BSIM3cgsb + BSIM3cbsb);
                BSIM3cqbb = -(BSIM3cqgb + BSIM3cqdb
                                    + BSIM3cqsb);

                gtau_drift = Math.Abs(param.BSIM3tconst * qcheq) * lp.ScalingFactor;
                T0 = param.BSIM3leffCV * param.BSIM3leffCV;
                gtau_diff = 16.0 * param.BSIM3u0temp * model.BSIM3vtm / T0
                            * lp.ScalingFactor;
                BSIM3gtau = gtau_drift + gtau_diff;
                if (BSIM3acnqsMod != 0)
                    BSIM3taunet = lp.ScalingFactor / BSIM3gtau;

            }

            if (model.BSIM3capMod == 0) /* code merge -JX */
            {
                cgdo = param.BSIM3cgdo;
                qgdo = param.BSIM3cgdo * vgd;
                cgso = param.BSIM3cgso;
                qgso = param.BSIM3cgso * vgs;
            }
            else if (model.BSIM3capMod == 1)
            {
                if (vgd < 0.0)
                {
                    T1 = Math.Sqrt(1.0 - 4.0 * vgd / param.BSIM3ckappa);
                    cgdo = param.BSIM3cgdo + param.BSIM3weffCV
                           * param.BSIM3cgdl / T1;
                    qgdo = param.BSIM3cgdo * vgd - param.BSIM3weffCV * 0.5
                           * param.BSIM3cgdl * param.BSIM3ckappa * (T1 - 1.0);
                }
                else
                {
                    cgdo = param.BSIM3cgdo + param.BSIM3weffCV
                           * param.BSIM3cgdl;
                    qgdo = (param.BSIM3weffCV * param.BSIM3cgdl
                            + param.BSIM3cgdo) * vgd;
                }

                if (vgs < 0.0)
                {
                    T1 = Math.Sqrt(1.0 - 4.0 * vgs / param.BSIM3ckappa);
                    cgso = param.BSIM3cgso + param.BSIM3weffCV
                           * param.BSIM3cgsl / T1;
                    qgso = param.BSIM3cgso * vgs - param.BSIM3weffCV * 0.5
                           * param.BSIM3cgsl * param.BSIM3ckappa * (T1 - 1.0);
                }
                else
                {
                    cgso = param.BSIM3cgso + param.BSIM3weffCV
                           * param.BSIM3cgsl;
                    qgso = (param.BSIM3weffCV * param.BSIM3cgsl
                            + param.BSIM3cgso) * vgs;
                }
            }
            else
            {
                T0 = vgd + DELTA_1;
                T1 = Math.Sqrt(T0 * T0 + 4.0 * DELTA_1);
                T2 = 0.5 * (T0 - T1);

                T3 = param.BSIM3weffCV * param.BSIM3cgdl;
                T4 = Math.Sqrt(1.0 - 4.0 * T2 / param.BSIM3ckappa);
                cgdo = param.BSIM3cgdo + T3 - T3 * (1.0 - 1.0 / T4)
                       * (0.5 - 0.5 * T0 / T1);
                qgdo = (param.BSIM3cgdo + T3) * vgd - T3 * (T2
                        + 0.5 * param.BSIM3ckappa * (T4 - 1.0));

                T0 = vgs + DELTA_1;
                T1 = Math.Sqrt(T0 * T0 + 4.0 * DELTA_1);
                T2 = 0.5 * (T0 - T1);
                T3 = param.BSIM3weffCV * param.BSIM3cgsl;
                T4 = Math.Sqrt(1.0 - 4.0 * T2 / param.BSIM3ckappa);
                cgso = param.BSIM3cgso + T3 - T3 * (1.0 - 1.0 / T4)
                       * (0.5 - 0.5 * T0 / T1);
                qgso = (param.BSIM3cgso + T3) * vgs - T3 * (T2
                        + 0.5 * param.BSIM3ckappa * (T4 - 1.0));
            }

            BSIM3cgdo = cgdo;
            BSIM3cgso = cgso;

            ag0 = integ.Deriv;
            if (BSIM3mode > 0)
            {
                if (BSIM3nqsMod == 0)
                {
                    gcggb = (BSIM3cggb + cgdo + cgso
                             + param.BSIM3cgbo) * ag0;
                    gcgdb = (BSIM3cgdb - cgdo) * ag0;
                    gcgsb = (BSIM3cgsb - cgso) * ag0;

                    gcdgb = (BSIM3cdgb - cgdo) * ag0;
                    gcddb = (BSIM3cddb + BSIM3capbd + cgdo) * ag0;
                    gcdsb = BSIM3cdsb * ag0;

                    gcsgb = -(BSIM3cggb + BSIM3cbgb
                              + BSIM3cdgb + cgso) * ag0;
                    gcsdb = -(BSIM3cgdb + BSIM3cbdb
                              + BSIM3cddb) * ag0;
                    gcssb = (BSIM3capbs + cgso - (BSIM3cgsb
                                                        + BSIM3cbsb + BSIM3cdsb)) * ag0;

                    gcbgb = (BSIM3cbgb - param.BSIM3cgbo) * ag0;
                    gcbdb = (BSIM3cbdb - BSIM3capbd) * ag0;
                    gcbsb = (BSIM3cbsb - BSIM3capbs) * ag0;

                    qgd = qgdo;
                    qgs = qgso;
                    qgb = param.BSIM3cgbo * vgb;
                    qgate += qgd + qgs + qgb;
                    qbulk -= qgb;
                    qdrn -= qgd;
                    qsrc = -(qgate + qbulk + qdrn);

                    ggtg = ggtd = ggtb = ggts = 0.0;
                    sxpart = 0.6;
                    dxpart = 0.4;
                    ddxpart_dVd = ddxpart_dVg = ddxpart_dVb = ddxpart_dVs = 0.0;
                    dsxpart_dVd = dsxpart_dVg = dsxpart_dVb = dsxpart_dVs = 0.0;
                }
                else
                {
                    if (qcheq > 0.0)
                        T0 = param.BSIM3tconst * qdef * lp.ScalingFactor;
                    else
                        T0 = -param.BSIM3tconst * qdef * lp.ScalingFactor;
                    ggtg = BSIM3gtg = T0 * BSIM3cqgb;
                    ggtd = BSIM3gtd = T0 * BSIM3cqdb;
                    ggts = BSIM3gts = T0 * BSIM3cqsb;
                    ggtb = BSIM3gtb = T0 * BSIM3cqbb;
                    gqdef = lp.ScalingFactor * ag0;

                    gcqgb = BSIM3cqgb * ag0;
                    gcqdb = BSIM3cqdb * ag0;
                    gcqsb = BSIM3cqsb * ag0;
                    gcqbb = BSIM3cqbb * ag0;

                    gcggb = (cgdo + cgso + param.BSIM3cgbo) * ag0;
                    gcgdb = -cgdo * ag0;
                    gcgsb = -cgso * ag0;

                    gcdgb = -cgdo * ag0;
                    gcddb = (BSIM3capbd + cgdo) * ag0;
                    gcdsb = 0.0;

                    gcsgb = -cgso * ag0;
                    gcsdb = 0.0;
                    gcssb = (BSIM3capbs + cgso) * ag0;

                    gcbgb = -param.BSIM3cgbo * ag0;
                    gcbdb = -BSIM3capbd * ag0;
                    gcbsb = -BSIM3capbs * ag0;

                    CoxWL = model.BSIM3cox * param.BSIM3weffCV
                            * param.BSIM3leffCV;
                    if (Math.Abs(qcheq) <= 1.0e-5 * CoxWL)
                    {
                        if (model.BSIM3xpart < 0.5)
                        {
                            dxpart = 0.4;
                        }
                        else if (model.BSIM3xpart > 0.5)
                        {
                            dxpart = 0.0;
                        }
                        else
                        {
                            dxpart = 0.5;
                        }
                        ddxpart_dVd = ddxpart_dVg = ddxpart_dVb
                                                    = ddxpart_dVs = 0.0;
                    }
                    else
                    {
                        dxpart = qdrn / qcheq;
                        Cdd = BSIM3cddb;
                        Csd = -(BSIM3cgdb + BSIM3cddb
                                + BSIM3cbdb);
                        ddxpart_dVd = (Cdd - dxpart * (Cdd + Csd)) / qcheq;
                        Cdg = BSIM3cdgb;
                        Csg = -(BSIM3cggb + BSIM3cdgb
                                + BSIM3cbgb);
                        ddxpart_dVg = (Cdg - dxpart * (Cdg + Csg)) / qcheq;

                        Cds = BSIM3cdsb;
                        Css = -(BSIM3cgsb + BSIM3cdsb
                                + BSIM3cbsb);
                        ddxpart_dVs = (Cds - dxpart * (Cds + Css)) / qcheq;

                        ddxpart_dVb = -(ddxpart_dVd + ddxpart_dVg + ddxpart_dVs);
                    }
                    sxpart = 1.0 - dxpart;
                    dsxpart_dVd = -ddxpart_dVd;
                    dsxpart_dVg = -ddxpart_dVg;
                    dsxpart_dVs = -ddxpart_dVs;
                    dsxpart_dVb = -(dsxpart_dVd + dsxpart_dVg + dsxpart_dVs);

                    qgd = qgdo;
                    qgs = qgso;
                    qgb = param.BSIM3cgbo * vgb;
                    qgate = qgd + qgs + qgb;
                    qbulk = -qgb;
                    qdrn = -qgd;
                    qsrc = -(qgate + qbulk + qdrn);
                }
            }
            else
            {
                if (BSIM3nqsMod == 0)
                {
                    gcggb = (BSIM3cggb + cgdo + cgso
                             + param.BSIM3cgbo) * ag0;
                    gcgdb = (BSIM3cgsb - cgdo) * ag0;
                    gcgsb = (BSIM3cgdb - cgso) * ag0;

                    gcdgb = -(BSIM3cggb + BSIM3cbgb
                              + BSIM3cdgb + cgdo) * ag0;
                    gcddb = (BSIM3capbd + cgdo - (BSIM3cgsb
                                                        + BSIM3cbsb + BSIM3cdsb)) * ag0;
                    gcdsb = -(BSIM3cgdb + BSIM3cbdb
                              + BSIM3cddb) * ag0;

                    gcsgb = (BSIM3cdgb - cgso) * ag0;
                    gcsdb = BSIM3cdsb * ag0;
                    gcssb = (BSIM3cddb + BSIM3capbs + cgso) * ag0;

                    gcbgb = (BSIM3cbgb - param.BSIM3cgbo) * ag0;
                    gcbdb = (BSIM3cbsb - BSIM3capbd) * ag0;
                    gcbsb = (BSIM3cbdb - BSIM3capbs) * ag0;

                    qgd = qgdo;
                    qgs = qgso;
                    qgb = param.BSIM3cgbo * vgb;
                    qgate += qgd + qgs + qgb;
                    qbulk -= qgb;
                    qsrc = qdrn - qgs;
                    qdrn = -(qgate + qbulk + qsrc);

                    ggtg = ggtd = ggtb = ggts = 0.0;
                    sxpart = 0.4;
                    dxpart = 0.6;
                    ddxpart_dVd = ddxpart_dVg = ddxpart_dVb = ddxpart_dVs = 0.0;
                    dsxpart_dVd = dsxpart_dVg = dsxpart_dVb = dsxpart_dVs = 0.0;
                }
                else
                {
                    if (qcheq > 0.0)
                        T0 = param.BSIM3tconst * qdef * lp.ScalingFactor;
                    else
                        T0 = -param.BSIM3tconst * qdef * lp.ScalingFactor;
                    ggtg = BSIM3gtg = T0 * BSIM3cqgb;
                    ggts = BSIM3gtd = T0 * BSIM3cqdb;
                    ggtd = BSIM3gts = T0 * BSIM3cqsb;
                    ggtb = BSIM3gtb = T0 * BSIM3cqbb;
                    gqdef = lp.ScalingFactor * ag0;

                    gcqgb = BSIM3cqgb * ag0;
                    gcqdb = BSIM3cqsb * ag0;
                    gcqsb = BSIM3cqdb * ag0;
                    gcqbb = BSIM3cqbb * ag0;

                    gcggb = (cgdo + cgso + param.BSIM3cgbo) * ag0;
                    gcgdb = -cgdo * ag0;
                    gcgsb = -cgso * ag0;

                    gcdgb = -cgdo * ag0;
                    gcddb = (BSIM3capbd + cgdo) * ag0;
                    gcdsb = 0.0;

                    gcsgb = -cgso * ag0;
                    gcsdb = 0.0;
                    gcssb = (BSIM3capbs + cgso) * ag0;

                    gcbgb = -param.BSIM3cgbo * ag0;
                    gcbdb = -BSIM3capbd * ag0;
                    gcbsb = -BSIM3capbs * ag0;

                    CoxWL = model.BSIM3cox * param.BSIM3weffCV
                            * param.BSIM3leffCV;
                    if (Math.Abs(qcheq) <= 1.0e-5 * CoxWL)
                    {
                        if (model.BSIM3xpart < 0.5)
                        {
                            sxpart = 0.4;
                        }
                        else if (model.BSIM3xpart > 0.5)
                        {
                            sxpart = 0.0;
                        }
                        else
                        {
                            sxpart = 0.5;
                        }
                        dsxpart_dVd = dsxpart_dVg = dsxpart_dVb
                                                    = dsxpart_dVs = 0.0;
                    }
                    else
                    {
                        sxpart = qdrn / qcheq;
                        Css = BSIM3cddb;
                        Cds = -(BSIM3cgdb + BSIM3cddb
                                + BSIM3cbdb);
                        dsxpart_dVs = (Css - sxpart * (Css + Cds)) / qcheq;
                        Csg = BSIM3cdgb;
                        Cdg = -(BSIM3cggb + BSIM3cdgb
                                + BSIM3cbgb);
                        dsxpart_dVg = (Csg - sxpart * (Csg + Cdg)) / qcheq;

                        Csd = BSIM3cdsb;
                        Cdd = -(BSIM3cgsb + BSIM3cdsb
                                + BSIM3cbsb);
                        dsxpart_dVd = (Csd - sxpart * (Csd + Cdd)) / qcheq;

                        dsxpart_dVb = -(dsxpart_dVd + dsxpart_dVg + dsxpart_dVs);
                    }
                    dxpart = 1.0 - sxpart;
                    ddxpart_dVd = -dsxpart_dVd;
                    ddxpart_dVg = -dsxpart_dVg;
                    ddxpart_dVs = -dsxpart_dVs;
                    ddxpart_dVb = -(ddxpart_dVd + ddxpart_dVg + ddxpart_dVs);

                    qgd = qgdo;
                    qgs = qgso;
                    qgb = param.BSIM3cgbo * vgb;
                    qgate = qgd + qgs + qgb;
                    qbulk = -qgb;
                    qsrc = -qgs;
                    qdrn = -(qgate + qbulk + qsrc);
                }
            }

            cqdef = cqcheq = 0.0;
            if (ByPass) goto line860;

            state.States[0][BSIM3state + BSIM3qg] = qgate;
            state.States[0][BSIM3state + BSIM3qd] = qdrn
                                                - state.States[0][BSIM3state + BSIM3qbd];
            state.States[0][BSIM3state + BSIM3qb] = qbulk
                                                + state.States[0][BSIM3state + BSIM3qbd]
                                                + state.States[0][BSIM3state + BSIM3qbs];

            if (BSIM3nqsMod != 0)
            {
                state.States[0][BSIM3state + BSIM3qcdump] = qdef * lp.ScalingFactor;
                state.States[0][BSIM3state + BSIM3qcheq] = qcheq;
            }

            /* store small signal parameters */
            if (ckt.Mode.HasFlag(Circuit.Modes.InitSmSig))
            {
                goto line1000;
            }
            if (!lp.ChargeComputationNeeded)
                goto line850;

            if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
            {
                state.States[1][BSIM3state + BSIM3qb] =
                    state.States[0][BSIM3state + BSIM3qb];
                state.States[1][BSIM3state + BSIM3qg] =
                    state.States[0][BSIM3state + BSIM3qg];
                state.States[1][BSIM3state + BSIM3qd] =
                    state.States[0][BSIM3state + BSIM3qd];
                if (BSIM3nqsMod != 0)
                {
                    state.States[1][BSIM3state + BSIM3qcheq] =
                        state.States[0][BSIM3state + BSIM3qcheq];
                    state.States[1][BSIM3state + BSIM3qcdump] =
                        state.States[0][BSIM3state + BSIM3qcdump];
                }
            }

            integ.Integrate(state, BSIM3state + BSIM3qb, 0.0);
            integ.Integrate(state, BSIM3state + BSIM3qg, 0.0);
            integ.Integrate(state, BSIM3state + BSIM3qd, 0.0);
            if (BSIM3nqsMod != 0)
            {
                integ.Integrate(state, BSIM3state + BSIM3qcdump, 0.0);
                integ.Integrate(state, BSIM3state + BSIM3qcheq, 0.0);
            }

            goto line860;

        line850:
            /* initialize to zero charge conductance and current */
            ceqqg = ceqqb = ceqqd = 0.0;
            cqcheq = cqdef = 0.0;

            gcdgb = gcddb = gcdsb = 0.0;
            gcsgb = gcsdb = gcssb = 0.0;
            gcggb = gcgdb = gcgsb = 0.0;
            gcbgb = gcbdb = gcbsb = 0.0;

            gqdef = gcqgb = gcqdb = gcqsb = gcqbb = 0.0;
            ggtg = ggtd = ggtb = ggts = 0.0;
            sxpart = (1.0 - (dxpart = (BSIM3mode > 0) ? 0.4 : 0.6));
            ddxpart_dVd = ddxpart_dVg = ddxpart_dVb = ddxpart_dVs = 0.0;
            dsxpart_dVd = dsxpart_dVg = dsxpart_dVb = dsxpart_dVs = 0.0;

            if (BSIM3nqsMod != 0)
                BSIM3gtau = 16.0 * param.BSIM3u0temp * model.BSIM3vtm
                                  / param.BSIM3leffCV / param.BSIM3leffCV
                                  * lp.ScalingFactor;
            else
                BSIM3gtau = 0.0;

            goto line900;

        line860:
            /* evaluate equivalent charge current */

            cqgate = state.States[0][BSIM3state + BSIM3cqg];
            cqbulk = state.States[0][BSIM3state + BSIM3cqb];
            cqdrn = state.States[0][BSIM3state + BSIM3cqd];

            ceqqg = cqgate - gcggb * vgb + gcgdb * vbd + gcgsb * vbs;
            ceqqb = cqbulk - gcbgb * vgb + gcbdb * vbd + gcbsb * vbs;
            ceqqd = cqdrn - gcdgb * vgb + gcddb * vbd + gcdsb * vbs;

            if (BSIM3nqsMod != 0)
            {
                T0 = ggtg * vgb - ggtd * vbd - ggts * vbs;
                ceqqg += T0;
                T1 = qdef * BSIM3gtau;
                ceqqd -= dxpart * T0 + T1 * (ddxpart_dVg * vgb - ddxpart_dVd
                                             * vbd - ddxpart_dVs * vbs);
                cqdef = state.States[0][BSIM3state + BSIM3cqcdump] - gqdef * qdef;
                cqcheq = state.States[0][BSIM3state + BSIM3cqcheq]
                         - (gcqgb * vgb - gcqdb * vbd - gcqsb * vbs) + T0;
            }

            if (ckt.Mode.HasFlag(Circuit.Modes.InitTran))
            {
                state.States[1][BSIM3state + BSIM3cqb] =
                    state.States[0][BSIM3state + BSIM3cqb];
                state.States[1][BSIM3state + BSIM3cqg] =
                    state.States[0][BSIM3state + BSIM3cqg];
                state.States[1][BSIM3state + BSIM3cqd] =
                    state.States[0][BSIM3state + BSIM3cqd];

                if (BSIM3nqsMod != 0)
                {
                    state.States[1][BSIM3state + BSIM3cqcheq] =
                        state.States[0][BSIM3state + BSIM3cqcheq];
                    state.States[1][BSIM3state + BSIM3cqcdump] =
                        state.States[0][BSIM3state + BSIM3cqcdump];
                }
            }

            /*
             *  load current vector
             */
        line900:

            if (BSIM3mode >= 0)
            {
                Gm = BSIM3gm;
                Gmbs = BSIM3gmbs;
                FwdSum = Gm + Gmbs;
                RevSum = 0.0;
                cdreq = model.BSIM3type * (cdrain - BSIM3gds * vds
                                            - Gm * vgs - Gmbs * vbs);

                ceqbd = -model.BSIM3type * (BSIM3csub
                                             - BSIM3gbds * vds - BSIM3gbgs * vgs
                                             - BSIM3gbbs * vbs);
                ceqbs = 0.0;

                gbbdp = -BSIM3gbds;
                gbbsp = (BSIM3gbds + BSIM3gbgs + BSIM3gbbs);

                gbdpg = BSIM3gbgs;
                gbdpdp = BSIM3gbds;
                gbdpb = BSIM3gbbs;
                gbdpsp = -(gbdpg + gbdpdp + gbdpb);

                gbspg = 0.0;
                gbspdp = 0.0;
                gbspb = 0.0;
                gbspsp = 0.0;
            }
            else
            {
                Gm = -BSIM3gm;
                Gmbs = -BSIM3gmbs;
                FwdSum = 0.0;
                RevSum = -(Gm + Gmbs);
                cdreq = -model.BSIM3type * (cdrain + BSIM3gds * vds
                                             + Gm * vgd + Gmbs * vbd);

                ceqbs = -model.BSIM3type * (BSIM3csub
                                             + BSIM3gbds * vds - BSIM3gbgs * vgd
                                             - BSIM3gbbs * vbd);
                ceqbd = 0.0;

                gbbsp = -BSIM3gbds;
                gbbdp = (BSIM3gbds + BSIM3gbgs + BSIM3gbbs);

                gbdpg = 0.0;
                gbdpsp = 0.0;
                gbdpb = 0.0;
                gbdpdp = 0.0;

                gbspg = BSIM3gbgs;
                gbspsp = BSIM3gbds;
                gbspb = BSIM3gbbs;
                gbspdp = -(gbspg + gbspsp + gbspb);
            }

            if (model.BSIM3type > 0)
            {
                ceqbs += (BSIM3cbs - BSIM3gbs * vbs);
                ceqbd += (BSIM3cbd - BSIM3gbd * vbd);
                /*
                    ceqqg = ceqqg;
                    ceqqb = ceqqb;
                    ceqqd = ceqqd;
                    cqdef = cqdef;
                    cqcheq = cqcheq;
                */
            }
            else
            {
                ceqbs -= (BSIM3cbs - BSIM3gbs * vbs);
                ceqbd -= (BSIM3cbd - BSIM3gbd * vbd);
                ceqqg = -ceqqg;
                ceqqb = -ceqqb;
                ceqqd = -ceqqd;
                cqdef = -cqdef;
                cqcheq = -cqcheq;
            }

            state.Rhs[BSIM3gNode] -= ceqqg;
            state.Rhs[BSIM3bNode] -= (ceqbs + ceqbd + ceqqb);
            state.Rhs[BSIM3dNodePrime] += (ceqbd - cdreq - ceqqd);
            state.Rhs[BSIM3sNodePrime] += (cdreq + ceqbs + ceqqg
                    + ceqqb + ceqqd);
            if (BSIM3nqsMod != 0)
                state.Rhs[BSIM3qNode] += (cqcheq - cqdef);

            /*
             *  load y matrix
             */

            T1 = qdef * BSIM3gtau;
            state.Matrix[BSIM3dNode, BSIM3dNode] += BSIM3drainConductance;
            state.Matrix[BSIM3gNode, BSIM3gNode] += gcggb - ggtg;
            state.Matrix[BSIM3sNode, BSIM3sNode] += BSIM3sourceConductance;
            state.Matrix[BSIM3bNode, BSIM3bNode] += BSIM3gbd + BSIM3gbs
                                    - gcbgb - gcbdb - gcbsb - BSIM3gbbs;
            state.Matrix[BSIM3dNodePrime, BSIM3dNodePrime] += BSIM3drainConductance
                                      + BSIM3gds + BSIM3gbd
                                      + RevSum + gcddb + dxpart * ggtd
                                      + T1 * ddxpart_dVd + gbdpdp;
            state.Matrix[BSIM3sNodePrime, BSIM3sNodePrime] += BSIM3sourceConductance
                                      + BSIM3gds + BSIM3gbs
                                      + FwdSum + gcssb + sxpart * ggts
                                      + T1 * dsxpart_dVs + gbspsp;
            state.Matrix[BSIM3dNode, BSIM3dNodePrime] -= BSIM3drainConductance;
            state.Matrix[BSIM3gNode, BSIM3bNode] -= gcggb + gcgdb + gcgsb + ggtb;
            state.Matrix[BSIM3gNode, BSIM3dNodePrime] += gcgdb - ggtd;
            state.Matrix[BSIM3gNode, BSIM3sNodePrime] += gcgsb - ggts;
            state.Matrix[BSIM3sNode, BSIM3sNodePrime] -= BSIM3sourceConductance;
            state.Matrix[BSIM3bNode, BSIM3gNode] += gcbgb - BSIM3gbgs;
            state.Matrix[BSIM3bNode, BSIM3dNodePrime] += gcbdb - BSIM3gbd + gbbdp;
            state.Matrix[BSIM3bNode, BSIM3sNodePrime] += gcbsb - BSIM3gbs + gbbsp;
            state.Matrix[BSIM3dNodePrime, BSIM3dNode] -= BSIM3drainConductance;
            state.Matrix[BSIM3dNodePrime, BSIM3gNode] += Gm + gcdgb + dxpart * ggtg
                                     + T1 * ddxpart_dVg + gbdpg;
            state.Matrix[BSIM3dNodePrime, BSIM3bNode] -= BSIM3gbd - Gmbs + gcdgb + gcddb
                                     + gcdsb - dxpart * ggtb
                                     - T1 * ddxpart_dVb - gbdpb;
            state.Matrix[BSIM3dNodePrime, BSIM3sNodePrime] -= BSIM3gds + FwdSum - gcdsb
                                      - dxpart * ggts - T1 * ddxpart_dVs - gbdpsp;
            state.Matrix[BSIM3sNodePrime, BSIM3gNode] += gcsgb - Gm + sxpart * ggtg
                                     + T1 * dsxpart_dVg + gbspg;
            state.Matrix[BSIM3sNodePrime, BSIM3sNode] -= BSIM3sourceConductance;
            state.Matrix[BSIM3sNodePrime, BSIM3bNode] -= BSIM3gbs + Gmbs + gcsgb + gcsdb
                                     + gcssb - sxpart * ggtb
                                     - T1 * dsxpart_dVb - gbspb;
            state.Matrix[BSIM3sNodePrime, BSIM3dNodePrime] -= BSIM3gds + RevSum - gcsdb
                                      - sxpart * ggtd - T1 * dsxpart_dVd - gbspdp;

            if (BSIM3nqsMod != 0)
            {
                state.Matrix[BSIM3qNode, BSIM3qNode] += (gqdef + BSIM3gtau);

                state.Matrix[BSIM3dNodePrime, BSIM3qNode] += (dxpart * BSIM3gtau);
                state.Matrix[BSIM3sNodePrime, BSIM3qNode] += (sxpart * BSIM3gtau);
                state.Matrix[BSIM3gNode, BSIM3qNode] -= BSIM3gtau;

                state.Matrix[BSIM3qNode, BSIM3gNode] += (ggtg - gcqgb);
                state.Matrix[BSIM3qNode, BSIM3dNodePrime] += (ggtd - gcqdb);
                state.Matrix[BSIM3qNode, BSIM3sNodePrime] += (ggts - gcqsb);
                state.Matrix[BSIM3qNode, BSIM3bNode] += (ggtb - gcqbb);
            }

        line1000:
            ;
        }

        /// <summary>
        /// Load the BSIM3 model in the circuit for AC simulation
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void AcLoad(Circuit ckt, Model mod, object p = null)
        {
            BSIM3Model model = mod as BSIM3Model;
            double xcggb, xcgdb, xcgsb, xcgbb, xcbgb, xcbdb, xcbsb, xcbbb;
            double xcdgb, xcddb, xcdsb, xcdbb, xcsgb, xcsdb, xcssb, xcsbb;
            double gdpr, gspr, gds, gbd, gbs, capbd, capbs, FwdSum, RevSum, Gm, Gmbs;
            double cggb, cgdb, cgsb, cbgb, cbdb, cbsb, cddb, cdgb, cdsb;
            double GSoverlapCap, GDoverlapCap, GBoverlapCap;
            double dxpart, sxpart, xgtg, xgtd, xgts, xgtb, xcqgb = 0.0, xcqdb = 0.0, xcqsb = 0.0, xcqbb = 0.0;
            double gbspsp, gbbdp, gbbsp, gbspg, gbspb;
            double gbspdp, gbdpdp, gbdpg, gbdpb, gbdpsp;
            double ddxpart_dVd, ddxpart_dVg, ddxpart_dVb, ddxpart_dVs;
            double dsxpart_dVd, dsxpart_dVg, dsxpart_dVb, dsxpart_dVs;
            double T1, CoxWL, qcheq, Cdg, Cdd, Cds, Csg, Csd, Css; // Cdb, Csb
            double ScalingFactor = 1.0e-9;

            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            if (BSIM3mode >= 0)
            {
                Gm = BSIM3gm;
                Gmbs = BSIM3gmbs;
                FwdSum = Gm + Gmbs;
                RevSum = 0.0;

                gbbdp = -BSIM3gbds;
                gbbsp = BSIM3gbds + BSIM3gbgs + BSIM3gbbs;

                gbdpg = BSIM3gbgs;
                gbdpdp = BSIM3gbds;
                gbdpb = BSIM3gbbs;
                gbdpsp = -(gbdpg + gbdpdp + gbdpb);

                gbspg = 0.0;
                gbspdp = 0.0;
                gbspb = 0.0;
                gbspsp = 0.0;

                if (BSIM3nqsMod == 0 && BSIM3acnqsMod == 0)
                {
                    cggb = BSIM3cggb;
                    cgsb = BSIM3cgsb;
                    cgdb = BSIM3cgdb;

                    cbgb = BSIM3cbgb;
                    cbsb = BSIM3cbsb;
                    cbdb = BSIM3cbdb;

                    cdgb = BSIM3cdgb;
                    cdsb = BSIM3cdsb;
                    cddb = BSIM3cddb;

                    xgtg = xgtd = xgts = xgtb = 0.0;
                    sxpart = 0.6;
                    dxpart = 0.4;
                    ddxpart_dVd = ddxpart_dVg = ddxpart_dVb
                        = ddxpart_dVs = 0.0;
                    dsxpart_dVd = dsxpart_dVg = dsxpart_dVb
                        = dsxpart_dVs = 0.0;
                }
                else
                {
                    cggb = cgdb = cgsb = 0.0;
                    cbgb = cbdb = cbsb = 0.0;
                    cdgb = cddb = cdsb = 0.0;

                    xgtg = BSIM3gtg;
                    xgtd = BSIM3gtd;
                    xgts = BSIM3gts;
                    xgtb = BSIM3gtb;

                    xcqgb = BSIM3cqgb;
                    xcqdb = BSIM3cqdb;
                    xcqsb = BSIM3cqsb;
                    xcqbb = BSIM3cqbb;

                    CoxWL = model.BSIM3cox * param.BSIM3weffCV
                                  * param.BSIM3leffCV;
                    qcheq = -(BSIM3qgate + BSIM3qbulk);
                    if (Math.Abs(qcheq) <= 1.0e-5 * CoxWL)
                    {
                        if (model.BSIM3xpart < 0.5)
                        {
                            dxpart = 0.4;
                        }
                        else if (model.BSIM3xpart > 0.5)
                        {
                            dxpart = 0.0;
                        }
                        else
                        {
                            dxpart = 0.5;
                        }
                        ddxpart_dVd = ddxpart_dVg = ddxpart_dVb
                            = ddxpart_dVs = 0.0;
                    }
                    else
                    {
                        dxpart = BSIM3qdrn / qcheq;
                        Cdd = BSIM3cddb;
                        Csd = -(BSIM3cgdb + BSIM3cddb
                        + BSIM3cbdb);
                        ddxpart_dVd = (Cdd - dxpart * (Cdd + Csd)) / qcheq;
                        Cdg = BSIM3cdgb;
                        Csg = -(BSIM3cggb + BSIM3cdgb
                        + BSIM3cbgb);
                        ddxpart_dVg = (Cdg - dxpart * (Cdg + Csg)) / qcheq;

                        Cds = BSIM3cdsb;
                        Css = -(BSIM3cgsb + BSIM3cdsb
                        + BSIM3cbsb);
                        ddxpart_dVs = (Cds - dxpart * (Cds + Css)) / qcheq;

                        ddxpart_dVb = -(ddxpart_dVd + ddxpart_dVg
                            + ddxpart_dVs);
                    }
                    sxpart = 1.0 - dxpart;
                    dsxpart_dVd = -ddxpart_dVd;
                    dsxpart_dVg = -ddxpart_dVg;
                    dsxpart_dVs = -ddxpart_dVs;
                    dsxpart_dVb = -(dsxpart_dVd + dsxpart_dVg + dsxpart_dVs);
                }
            }
            else
            {
                Gm = -BSIM3gm;
                Gmbs = -BSIM3gmbs;
                FwdSum = 0.0;
                RevSum = -(Gm + Gmbs);

                gbbsp = -BSIM3gbds;
                gbbdp = BSIM3gbds + BSIM3gbgs + BSIM3gbbs;

                gbdpg = 0.0;
                gbdpsp = 0.0;
                gbdpb = 0.0;
                gbdpdp = 0.0;

                gbspg = BSIM3gbgs;
                gbspsp = BSIM3gbds;
                gbspb = BSIM3gbbs;
                gbspdp = -(gbspg + gbspsp + gbspb);

                if (BSIM3nqsMod == 0 && BSIM3acnqsMod == 0)
                {
                    cggb = BSIM3cggb;
                    cgsb = BSIM3cgdb;
                    cgdb = BSIM3cgsb;

                    cbgb = BSIM3cbgb;
                    cbsb = BSIM3cbdb;
                    cbdb = BSIM3cbsb;

                    cdgb = -(BSIM3cdgb + cggb + cbgb);
                    cdsb = -(BSIM3cddb + cgsb + cbsb);
                    cddb = -(BSIM3cdsb + cgdb + cbdb);

                    xgtg = xgtd = xgts = xgtb = 0.0;
                    sxpart = 0.4;
                    dxpart = 0.6;
                    ddxpart_dVd = ddxpart_dVg = ddxpart_dVb
                        = ddxpart_dVs = 0.0;
                    dsxpart_dVd = dsxpart_dVg = dsxpart_dVb
                        = dsxpart_dVs = 0.0;
                }
                else
                {
                    cggb = cgdb = cgsb = 0.0;
                    cbgb = cbdb = cbsb = 0.0;
                    cdgb = cddb = cdsb = 0.0;

                    xgtg = BSIM3gtg;
                    xgtd = BSIM3gts;
                    xgts = BSIM3gtd;
                    xgtb = BSIM3gtb;

                    xcqgb = BSIM3cqgb;
                    xcqdb = BSIM3cqsb;
                    xcqsb = BSIM3cqdb;
                    xcqbb = BSIM3cqbb;

                    CoxWL = model.BSIM3cox * param.BSIM3weffCV
                                  * param.BSIM3leffCV;
                    qcheq = -(BSIM3qgate + BSIM3qbulk);
                    if (Math.Abs(qcheq) <= 1.0e-5 * CoxWL)
                    {
                        if (model.BSIM3xpart < 0.5)
                        {
                            sxpart = 0.4;
                        }
                        else if (model.BSIM3xpart > 0.5)
                        {
                            sxpart = 0.0;
                        }
                        else
                        {
                            sxpart = 0.5;
                        }
                        dsxpart_dVd = dsxpart_dVg = dsxpart_dVb
                            = dsxpart_dVs = 0.0;
                    }
                    else
                    {
                        sxpart = BSIM3qdrn / qcheq;
                        Css = BSIM3cddb;
                        Cds = -(BSIM3cgdb + BSIM3cddb
                        + BSIM3cbdb);
                        dsxpart_dVs = (Css - sxpart * (Css + Cds)) / qcheq;
                        Csg = BSIM3cdgb;
                        Cdg = -(BSIM3cggb + BSIM3cdgb
                        + BSIM3cbgb);
                        dsxpart_dVg = (Csg - sxpart * (Csg + Cdg)) / qcheq;

                        Csd = BSIM3cdsb;
                        Cdd = -(BSIM3cgsb + BSIM3cdsb
                        + BSIM3cbsb);
                        dsxpart_dVd = (Csd - sxpart * (Csd + Cdd)) / qcheq;

                        dsxpart_dVb = -(dsxpart_dVd + dsxpart_dVg
                            + dsxpart_dVs);
                    }
                    dxpart = 1.0 - sxpart;
                    ddxpart_dVd = -dsxpart_dVd;
                    ddxpart_dVg = -dsxpart_dVg;
                    ddxpart_dVs = -dsxpart_dVs;
                    ddxpart_dVb = -(ddxpart_dVd + ddxpart_dVg + ddxpart_dVs);
                }
            }


            T1 = cstate.States[0][BSIM3state + BSIM3qdef] * BSIM3gtau;
            gdpr = BSIM3drainConductance;
            gspr = BSIM3sourceConductance;
            gds = BSIM3gds;
            gbd = BSIM3gbd;
            gbs = BSIM3gbs;
            capbd = BSIM3capbd;
            capbs = BSIM3capbs;

            GSoverlapCap = BSIM3cgso;
            GDoverlapCap = BSIM3cgdo;
            GBoverlapCap = param.BSIM3cgbo;

            xcdgb = (cdgb - GDoverlapCap);
            xcddb = (cddb + capbd + GDoverlapCap);
            xcdsb = cdsb;
            xcdbb = -(xcdgb + xcddb + xcdsb);
            xcsgb = -(cggb + cbgb + cdgb + GSoverlapCap);
            xcsdb = -(cgdb + cbdb + cddb);
            xcssb = (capbs + GSoverlapCap - (cgsb + cbsb + cdsb));
            xcsbb = -(xcsgb + xcsdb + xcssb);
            xcggb = (cggb + GDoverlapCap + GSoverlapCap + GBoverlapCap);
            xcgdb = (cgdb - GDoverlapCap);
            xcgsb = (cgsb - GSoverlapCap);
            xcgbb = -(xcggb + xcgdb + xcgsb);
            xcbgb = (cbgb - GBoverlapCap);
            xcbdb = (cbdb - capbd);
            xcbsb = (cbsb - capbs);
            xcbbb = -(xcbgb + xcbdb + xcbsb);

            state.Matrix[BSIM3gNode, BSIM3gNode] += xcggb * state.Omega;
            state.Matrix[BSIM3bNode, BSIM3bNode] += xcbbb * state.Omega;
            state.Matrix[BSIM3dNodePrime, BSIM3dNodePrime] += xcddb * state.Omega;
            state.Matrix[BSIM3sNodePrime, BSIM3sNodePrime] += xcssb * state.Omega;

            state.Matrix[BSIM3gNode, BSIM3bNode] += xcgbb * state.Omega;
            state.Matrix[BSIM3gNode, BSIM3dNodePrime] += xcgdb * state.Omega;
            state.Matrix[BSIM3gNode, BSIM3sNodePrime] += xcgsb * state.Omega;

            state.Matrix[BSIM3bNode, BSIM3gNode] += xcbgb * state.Omega;
            state.Matrix[BSIM3bNode, BSIM3dNodePrime] += xcbdb * state.Omega;
            state.Matrix[BSIM3bNode, BSIM3sNodePrime] += xcbsb * state.Omega;

            state.Matrix[BSIM3dNodePrime, BSIM3gNode] += xcdgb * state.Omega;
            state.Matrix[BSIM3dNodePrime, BSIM3bNode] += xcdbb * state.Omega;
            state.Matrix[BSIM3dNodePrime, BSIM3sNodePrime] += xcdsb * state.Omega;

            state.Matrix[BSIM3sNodePrime, BSIM3gNode] += xcsgb * state.Omega;
            state.Matrix[BSIM3sNodePrime, BSIM3bNode] += xcsbb * state.Omega;
            state.Matrix[BSIM3sNodePrime, BSIM3dNodePrime] += xcsdb * state.Omega;

            state.Matrix[BSIM3dNode, BSIM3dNode] += gdpr;
            state.Matrix[BSIM3dNode, BSIM3dNodePrime] -= gdpr;
            state.Matrix[BSIM3dNodePrime, BSIM3dNode] -= gdpr;

            state.Matrix[BSIM3sNode, BSIM3sNode] += gspr;
            state.Matrix[BSIM3sNode, BSIM3sNodePrime] -= gspr;
            state.Matrix[BSIM3sNodePrime, BSIM3sNode] -= gspr;

            state.Matrix[BSIM3bNode, BSIM3gNode] -= BSIM3gbgs;
            state.Matrix[BSIM3bNode, BSIM3bNode] += gbd + gbs - BSIM3gbbs;
            state.Matrix[BSIM3bNode, BSIM3dNodePrime] -= gbd - gbbdp;
            state.Matrix[BSIM3bNode, BSIM3sNodePrime] -= gbs - gbbsp;

            state.Matrix[BSIM3dNodePrime, BSIM3gNode] += Gm + dxpart * xgtg + T1 * ddxpart_dVg + gbdpg;
            state.Matrix[BSIM3dNodePrime, BSIM3dNodePrime] += gdpr + gds + gbd + RevSum + dxpart * xgtd + T1 * ddxpart_dVd + gbdpdp;
            state.Matrix[BSIM3dNodePrime, BSIM3sNodePrime] -= gds + FwdSum - dxpart * xgts - T1 * ddxpart_dVs - gbdpsp;
            state.Matrix[BSIM3dNodePrime, BSIM3bNode] -= gbd - Gmbs - dxpart * xgtb - T1 * ddxpart_dVb - gbdpb;

            state.Matrix[BSIM3sNodePrime, BSIM3gNode] -= Gm - sxpart * xgtg - T1 * dsxpart_dVg - gbspg;
            state.Matrix[BSIM3sNodePrime, BSIM3sNodePrime] += gspr + gds + gbs + FwdSum + sxpart * xgts + T1 * dsxpart_dVs + gbspsp;
            state.Matrix[BSIM3sNodePrime, BSIM3bNode] -= gbs + Gmbs - sxpart * xgtb - T1 * dsxpart_dVb - gbspb;
            state.Matrix[BSIM3sNodePrime, BSIM3dNodePrime] -= gds + RevSum - sxpart * xgtd - T1 * dsxpart_dVd - gbspdp;

            state.Matrix[BSIM3gNode, BSIM3gNode] -= xgtg;
            state.Matrix[BSIM3gNode, BSIM3bNode] -= xgtb;
            state.Matrix[BSIM3gNode, BSIM3dNodePrime] -= xgtd;
            state.Matrix[BSIM3gNode, BSIM3sNodePrime] -= xgts;

            if (BSIM3nqsMod != 0 || BSIM3acnqsMod != 0)
            {
                state.Matrix[BSIM3qNode, BSIM3qNode] += ScalingFactor * state.Omega;
                state.Matrix[BSIM3qNode, BSIM3gNode] -= xcqgb * state.Omega;
                state.Matrix[BSIM3qNode, BSIM3dNodePrime] -= xcqdb * state.Omega;
                state.Matrix[BSIM3qNode, BSIM3bNode] -= xcqbb * state.Omega;
                state.Matrix[BSIM3qNode, BSIM3sNodePrime] -= xcqsb * state.Omega;

                state.Matrix[BSIM3gNode, BSIM3qNode] -= BSIM3gtau;
                state.Matrix[BSIM3dNodePrime, BSIM3qNode] += dxpart * BSIM3gtau;
                state.Matrix[BSIM3sNodePrime, BSIM3qNode] += sxpart * BSIM3gtau;

                state.Matrix[BSIM3qNode, BSIM3qNode] += BSIM3gtau;
                state.Matrix[BSIM3qNode, BSIM3gNode] += xgtg;
                state.Matrix[BSIM3qNode, BSIM3dNodePrime] += xgtd;
                state.Matrix[BSIM3qNode, BSIM3bNode] += xgtb;
                state.Matrix[BSIM3qNode, BSIM3sNodePrime] += xgts;
            }
        }
    }
}
